The physical origin of long gas depletion times in galaxies

TL;DR

This paper introduces a model explaining why galaxies have long gas depletion times, emphasizing the role of feedback and multiple gas cycling through star-forming states, without assuming equilibrium.

Contribution

The model provides a non-equilibrium framework for understanding long depletion times and explains how feedback self-regulates star formation rates in galaxies.

Findings

Gas evolves quickly from non-star-forming to star-forming states.

Star formation rate is limited by the fraction of gas converted before dispersal.

The model reproduces observed Kennicutt-Schmidt relations for molecular and atomic gas.

Abstract

We present a model that explains why galaxies form stars on a time scale significantly longer than the time scales of processes governing the evolution of interstellar gas. We show that gas evolves from a non-star-forming to a star-forming state on a relatively short time scale and thus the rate of this evolution does not limit the star formation rate. Instead, the star formation rate is limited because only a small fraction of star-forming gas is converted into stars before star-forming regions are dispersed by feedback and dynamical processes. Thus, gas cycles into and out of star-forming state multiple times, which results in a long time scale on which galaxies convert gas into stars. Our model does not rely on the assumption of equilibrium and can be used to interpret trends of depletion times with the properties of observed galaxies and the parameters of star formation and feedback recipes in simulations. In particular, the model explains how feedback self-regulates the star formation rate in simulations and makes it insensitive to the local star formation efficiency. We illustrate our model using the results of an isolated $L_*$-sized galaxy simulation that reproduces the observed Kennicutt-Schmidt relation for both molecular and atomic gas. Interestingly, the relation for molecular gas is almost linear on kiloparsec scales, although a nonlinear relation is adopted in simulation cells. We discuss how a linear relation emerges from non-self-similar scaling of the gas density PDF with the average gas surface density.

Equations

$${\tau }_{\mathrm{dep}}\equiv \frac{M_{\mathrm{g}}}{{\dot{M}}_{\star }}=\frac{{\tau }_{\star }}{f_{\mathrm{sf}}},$$

$$t_{\mathrm{dep}}=N_{\mathrm{c}}(t_{\mathrm{nsf}}+t_{\mathrm{sf}})=N_{\mathrm{c}}{t}_{\mathrm{nsf}}+{\tau }_{\star }=\left(\frac{t_{\mathrm{nsf}}}{t_{\mathrm{sf}}}+1\right){\tau }_{\star }.$$

$${\tau }_{\mathrm{dep}}\equiv \frac{M_{\mathrm{g}}}{⟨{\dot{M}}_{\star }{⟩}_{\Delta t}}=\frac{M_{\mathrm{g}}}{⟨F_{\mathrm{sf}}{⟩}_{\Delta t}-⟨{\dot{M}}_{\mathrm{sf}}{⟩}_{\Delta t}}.$$

$$F_{\mathrm{sf}}=F_{+}-F_{-},$$

$$F_{+}$$

$$F_{-}$$

$$F_{-,\mathrm{fb}}\equiv \xi {\dot{M}}_{\star }=M_{\mathrm{g}}{f}_{\mathrm{sf}}\frac{\xi }{{\tau }_{\star }},$$

$$|{\dot{M}}_{\mathrm{sf}}|\equiv \frac{M_{\mathrm{sf}}}{{\tau }_{\mathrm{e},\mathrm{sf}}}.$$

$${\tau }_{\mathrm{dep}}=\left(1+\xi +\frac{{\tau }_{\star }}{{\tau }_{-,\mathrm{d}}}\pm \frac{{\tau }_{\star }}{{\tau }_{\mathrm{e},\mathrm{sf}}}\right){\tau }_{+}+{\tau }_{\star },$$

$$N_{\mathrm{c}}={\tau }_{\star }\left(\frac{1}{{\tau }_{\star }}+\frac{1}{{\tau }_{-}}\pm \frac{1}{{\tau }_{\mathrm{e},\mathrm{sf}}}\right)=1+\xi +\frac{{\tau }_{\star }}{{\tau }_{-,\mathrm{d}}}\pm \frac{{\tau }_{\star }}{{\tau }_{\mathrm{e},\mathrm{sf}}},$$

$${\alpha }_{\mathrm{vir}}\equiv \frac{5{\sigma }_{1D}^{2}R}{GM}\approx 9.35\frac{({\sigma }_{\mathrm{tot}}/10\mathrm{km}{\mathrm{s}}^{-1}{)}^2}{(n/100\mathrm{c}{\mathrm{m}}^{-3})(\Delta /40\mathrm{pc}{)}^2}.$$

$${\sigma }_{\mathrm{tot}}=\sqrt{{\sigma }_{\mathrm{t}}^2+c_{\mathrm{s}}^2}.$$

$${\tau }_{\alpha ,\gtrless 0}\equiv {⟨\left|\frac{d{\log }_{10}{\alpha }_{\mathrm{vir}}}{dt}\right|⟩}_{\gtrless 0}^{-1},$$

$${\tau }_{\mathrm{arm}}\sim \frac{2\pi R}{m{V}_{\mathrm{gas}}}\sim 80\mathrm{Myr}\frac{(R/8\mathrm{kpc})}{(m/6)(V_{\mathrm{gas}}/100\mathrm{km}{\mathrm{s}}^{-1})},$$

$${\tau }_{\mathrm{shell}}\sim \frac{L}{v_{\mathrm{shell}}}\sim 50\mathrm{Myr}\frac{(L/1\mathrm{kpc})}{(v_{\mathrm{shell}}/20\mathrm{km}{\mathrm{s}}^{-1})},$$

$${\tau }_{\mathrm{grav}}\sim \frac{v_z}{\pi G{\Sigma }_{\mathrm{tot}}}\sim 20\mathrm{Myr}\frac{(v_z/50\mathrm{km}{\mathrm{s}}^{-1})}{({\Sigma }_{\mathrm{tot}}/200{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-2})},$$

$$N_{\mathrm{c}}^{-1}=\int \frac{dt}{t_{\star }}=\int {ϵ}_{\mathrm{ff}}\frac{dt}{t_{\mathrm{ff}}},$$

$${\tau }_{\mathrm{dep},{\mathrm{H}}_2}\equiv \frac{{\Sigma }_{{\mathrm{H}}_2}}{{\dot{\Sigma }}_{\star }}=\frac{{\tau }_{\star }}{f_{\mathrm{sf},{\mathrm{H}}_2}},$$

Full text

††slugcomment: Accepted for publication in The Astrophysical Journal

1. Introduction

One of the widely recognized basic facts about observed galaxies is that they convert gas into stars inefficiently. Star formation rates of galaxies are surprisingly low, given the total mass and density of their interstellar gas and the expected time scales of processes driving star formation.

As an example, the star formation rate (SFR) of the Milky Way (MW) is $\dot{M}_{\star}\sim 1-2{\rm~{}M_{\odot}~{}yr^{-1}}$, while its gas mass is ${M}_{\rm g}\sim 10^{10}\ {\rm~{}M_{\odot}}$. Thus, the time scale at which the Galaxy would deplete its gas while forming stars at the current rate is $\tau_{\rm dep}\equiv{M}_{\rm g}/\dot{M}_{\star}\sim 5-10{\rm~{}Gyr}$. The depletion time scales of a population of normal star-forming galaxies are comparable and span a range of $\sim 2-10$ Gyr (Kennicutt.1989; Kennicutt.1998; Bigiel.etal.2008). The denser molecular phase of the interstellar medium (ISM) is depleted on a similarly long time scale of $\tau_{\rm dep,H_{2}}\equiv M_{\rm H_{2}}/\dot{M}_{\star}\sim 1-3$ Gyrs (Kennicutt.1989; Kennicutt.1998; wong_blitz02; Bigiel.etal.2008; leroy_etal08; Leroy.etal.2013).

Compared to the time scales of any dynamical processes that are potentially relevant for star formation, the observed gas depletion times are very long indeed. For example, the orbital period of gas at the solar radius is $t_{\rm orb}\sim 200$ Myr, and the MW is thus depleting its gas on the timescale of $\sim 25-50$ such periods. On average, galaxies deplete their gas on a timescale of $\sim 10-20$ orbital periods (Kennicutt.1998; wong_blitz02; leroy_etal08; daddi_etal10).

The orbital period, $t_{\rm orb}$, is probably the longest of the relevant dynamical timescales one can think of. For example, the turbulent crossing time is usually $t_{\rm cross}=h/\sigma\sim 10-30{\rm~{}Myr}$, where $\sigma\gtrsim 10{\rm~{}km~{}s^{-1}}$ is the velocity dispersion of gas in galactic disks and $h\sim 100-300{\rm~{}pc}$ is the disk scale height in the inner regions of galaxies. The free-fall time at the mean or midplane density, $\rho_{0}$, of galaxies spans a similar range: $t_{\rm ff,0}\equiv\sqrt{3\pi/32G\rho_{0}}\sim 10-50{\rm~{}Myr}$. The timescale of molecular cloud collisions is $\lesssim 20{\rm~{}Myr}$ (e.g., Tan.2000). A given gas mass encounters a spiral arm on a timescale of $t_{\rm arm}\sim 2\pi/(m[\Omega(R)-\Omega_{\rm p}])$, where $\Omega(R)=V_{\rm rot}/R$ is the angular frequency of gas rotation, $\Omega_{\rm p}$ is the pattern speed of spiral arms, and $m$ is the number of spiral arms. This timescale is $t_{\rm arm}\sim 50-200{\rm~{}Myr}$, if we assume $\Omega_{\rm p}\sim 20\rm{\rm~{}km~{}s^{-1}}\,kpc^{-1}$ (e.g., Bissantz.etal.2003), $m\sim 2-4$ (e.g., Davis.etal.2015) and $V_{\rm rot}\sim 220{\rm~{}km~{}s^{-1}}$ typically derived for MW-like galaxies. Numerical simulations of gaseous galactic disks show that star-forming molecular clouds may form on even shorter timescales of a few tens of Myrs (Dobbs.etal.2012; Dobbs.etal.2015).

In addition to being slow on global galactic scales, star formation is inefficient even in dense molecular star-forming regions, which convert only $\lesssim 1-10\%$ of gas into stars per local free-fall time (Zuckerman.Evans.1974; zuckerman_palmer74; Krumholz.Tan.2007; Krumholz.etal.2012; Evans.etal.2014; Lee.etal.2016; Heyer.et.al.2016). Such low efficiency arises because only $\sim 0.1-10\%$ of the dense gas is self-gravitating and collapsing into stars (Froebrich.Rowles.2010).

However, the inefficiency of star-forming regions alone cannot explain long global depletion times. The local depletion time in observed star-forming regions is $t_{\star}\sim 40-500{\rm~{}Myr}$ (e.g., Evans.etal.2009; Evans.etal.2014; Lada.etal.2010; Lada.etal.2012; Heiderman.etal.2010; Gutermuth.etal.2011; Schruba.etal.2017). Thus, although the scatter is significant, typical values of $t_{\star}$ are considerably smaller than the global depletion time of molecular gas, $\tau_{\rm dep,H_{2}}\sim 1-3{\rm~{}Gyr}$.

The large scatter in depletion times measured on small scales and the difference between local and global depletion time values indicate that only a fraction of molecular gas is actively forming stars at any given moment. Indeed, the global depletion time can be expressed as

[TABLE]

where ${M}_{\rm g}$ is the total gas mass of the galaxy; $\tau_{\star}\equiv{M}_{\rm sf}/\dot{M}_{\star}=\langle 1/t_{\star}\rangle_{\rm sf}^{-1}$ is the mass-weighted average over the depletion time distribution in star-forming regions, $t_{\star}$; and $f_{\rm sf}\equiv{M}_{\rm sf}/{M}_{\rm g}$ is the gas mass fraction in actively star-forming regions. A similar expression can be written for the global depletion time of molecular gas, $\tau_{\rm dep,H_{2}}$, via a corresponding star-forming fraction $f_{\rm sf,H_{2}}\equiv{M}_{\rm sf}/M_{\rm H_{2}}$.

Thus, the depletion time measured on larger scales is longer than that in star-forming regions because, as the scale increases, more of non-star-forming gas is incorporated in the gas mass estimate. Likewise, when depletion time is estimated on larger scales, the scatter in $\tau_{\rm dep}$ decreases as we average over the distribution of local $t_{\star}$. This is indeed observed (schruba_etal10; Schruba.etal.2017), although some of the obtained variation may be due to observational effects (Feldmann.etal.2011; Kruijssen.Longmore.2014).

The small scatter in the global depletion time in observed galaxies is manifested in the Kennicutt-Schmidt relation (KSR) between the surface density of gas and the star formation rate (Schmidt.1959; Kennicutt.1989; Kennicutt.1998, see also Sanduleak.1969; Madore.etal.1974). The relation is particularly tight and close to linear when only molecular hydrogen is used to estimate the surface density of gas (wong_blitz02; Bigiel.etal.2008; Leroy.etal.2013). Thus, a comprehensive model for the global depletion time must explain not only its value, i.e., the normalization of the KSR, but also both the scatter and the shape of the KSR on different scales.

A number of useful global star formation frameworks and models have been developed over the last three decades (e.g., Wyse.Silk.1989; Silk.1997; Tan.2000; Elmegreen.2002; Krumholz.McKee.2005; Li.etal.2005; Krumholz.Thompson.2007; Saitoh.etal.2008; Krumholz.etal.2009; Silk.Norman.2009; Ostriker.et.al.2010; Ostriker.Shetty.2011; Renaud.etal.2012; FaucherGiguere.etal.2013; Federrath.2013; Elmegreen.2015; Salim.etal.2015). Many of these models consider the physical processes shaping the form of the KSR while treating its normalization as a flexible constant. Other models also consider the physical origin of the normalization and long depletion timescale.

One class of the latter models associates long depletion times with the fraction of gas in dense, self-gravitating regions of cold, supersonic molecular clouds with the log-normal gas density PDF (Elmegreen.2002; Krumholz.McKee.2005; Krumholz.etal.2012). Such models, however, assume that all of the molecular gas is in “virialized” star-forming molecular clouds and that the star formation efficiency in these clouds sets the global depletion time. This assumption, which has also been frequently adopted in galaxy simulations (e.g., Robertson.Kravtsov.2008; Gnedin.etal.2009; Christensen.etal.2012; Kuhlen.etal.2012), is at odds with a growing number of observations indicating that the depletion time of star-forming molecular gas is in general considerably shorter than the global depletion time of all molecular gas, $\tau_{\rm dep,H_{2}}$, estimated on $\gtrsim$ kpc scales. Moreover, models and simulations of star formation in supersonic turbulent clouds show that the local efficiency of star formation is primarily a strong function of the virial parameter of the region, not just its density, temperature, and molecular fraction (e.g., Krumholz.McKee.2005; Padoan.etal.2012; Padoan.etal.2017), while the virial parameter can span a wide range of values (Dobbs.etal.2011; Semenov.etal.2016).

Some models derive the Kennicutt-Schmidt relation and its normalization by assuming that stellar feedback regulates ISM turbulence so as to maintain vertical and/or Toomre.1964 equilibrium within gaseous disks (Ostriker.Shetty.2011; FaucherGiguere.etal.2013; Hayward.Hopkins.2017). However, it is not clear a priori why equilibrium can generically be expected in galaxies as a whole or in kiloparsec-scale patches and why the star formation rate does not instead reach values at which gas is driven out in a wind. Moreover, it is still debated whether the turbulence within galactic disks is mainly driven by stellar feedback or by gravitational instabilities (e.g., Krumholz.Burkhart.2016).

Saitoh.etal.2008 argued that SFR is controlled by the rate at which gas is supplied from the general ISM to the star-forming state, which makes it insensitive to the local efficiency of star formation. However, these authors measured the timescale at which gas is supplied to the star-forming state to be $\sim 100$ Myr and did not explain how this timescale relates to the much longer observed depletion times of total gas, $\tau_{\rm dep}\sim 2-10$ Gyrs.

In this paper, we aim to clarify the origin of the observed long gas depletion timescale in galaxies, taking into account both the inefficiency of star formation in star-forming clouds and the fact that not all of the molecular gas is actively forming stars. To this end, in Section 2 we present a simple conceptual framework that views the ISM as a highly dynamic medium in which gas evolves between non-star-forming and star-forming states. The key aspect of our framework is that it considers global gas depletion as a result of gas evolution that is driven by processes with associated characteristic timescales. This approach is conceptually similar to the framework of Kruijssen.Longmore.2014, developed to explain and interpret the scatter and possible biases in observational measurements of the Kennicutt-Schmidt relation on different scales.

We illustrate our framework using a realistic simulation of a galactic disk that reproduces the observed depletion time and the Kennicutt-Schmidt relation. We describe the simulation in Section 3. In Section 4 we demonstrate that the long global depletion time originates from the rapid cycling of ISM gas between non-star-forming and star-forming states on timescales $\sim 20-100$ Myr, in accord with the above estimates and the results of previous galactic disk simulations. On each such cycle, only a small fraction of the gas mass is converted into stars and thus $\tau_{\rm dep}$ is long because a large number of such cycles would be required to deplete all available gas. We analyze the processes driving the rapid gas evolution and also use simulations to shed light on the reason why the depletion time of molecular gas is nearly independent of the gas surface density. We discuss our results and summarize our conclusions in Sections 5 and 6.

2. A model for gas depletion time

The interstellar gas in galaxies is a multiphase, dynamic medium spanning several orders of magnitude in density and temperature. To get a sense of the processes affecting the gas evolution in such a medium, we consider the evolution of individual gas parcels, massless tracers of gas flows in the ISM. One can think of a representative set of the ISM atoms as such tracers. At any given time, the local environment of such tracers can be estimated by averaging gas properties on some scale $l$ around the position of each tracer. The gas around tracers will expand and contract under the influence of dynamical processes, such as turbulence driven by gravitational instabilities and stellar feedback. Therefore, during evolution over a sufficiently long timescale, the conditions around each gas parcel can evolve between the states of long and short depletion time, $t_{\star}$, once or over many transition cycles. Such cycling of a gas parcel is schematically shown in Figure 1.

The probability density per unit time for a parcel to be converted into a star is given by $1/t_{\star}$ and one can define the depletion time for a single parcel as the time required for the integrated probability to reach unity. In what follows, for conceptual simplicity, we will adopt a sharp threshold, $t_{\star,\rm max}$, separating non-star-forming, $t_{\star}>t_{\star,\rm max}$, and actively star-forming, $t_{\star}<t_{\star,\rm max}$, gas states. For a given distribution of $t_{\star}$, this threshold can be chosen in such a way that regions with $t_{\star}<t_{\star,\rm max}$ include most of the total star formation.

The duration of a single cycle between the successive stages when the parcel’s environment is in the star-forming state equals to the sum of the time spent in the non-star-forming, $t_{\rm nsf}$, and star-forming, $t_{\rm sf}$, stages. If we denote the average depletion time of the parcel during the star-forming stage as $\tau_{\star}$, the parcel will have a probability of $t_{\rm sf}/\tau_{\star}$ to be converted into a star during the entire cycle. In other words, $N_{\rm c}=\tau_{\star}/t_{\rm sf}$ such cycles would be needed for the parcel to be incorporated into a star. Its depletion time can thus be written as

[TABLE]

Hence, the depletion time of a gas parcel is always longer than $\tau_{\star}$ and it can be long because star formation during $t_{\rm sf}$ is inherently inefficient, i.e., $\tau_{\star}$ is long, and/or because only a small fraction of the evolution cycle is spent in the star-forming state, i.e., $t_{\rm nsf}/t_{\rm sf}$ is large.

A group of parcels in a patch of the ISM has a distribution of $\tau_{\star}$, $t_{\rm sf}$, $t_{\rm nsf}$ and the collective depletion time of the group is the average $\tau_{\rm dep}=\langle 1/t_{\rm dep}\rangle^{-1}$ over these distributions. It is clear that if the depletion times of individual parcels are long, $\tau_{\rm dep}$ will also be long.

In practice, the $\tau_{\rm dep}$ of an ISM patch is estimated from the instantaneous mass of gas, ${M}_{\rm g}$, and young stars formed over a time interval $\Delta t$, ${M}_{\star}(<\Delta t)$, such that the average star formation rate is $\langle\dot{M}_{\star}\rangle_{\Delta t}={M}_{\star}(<\Delta t)/\Delta t$ and the depletion time is defined as $\tau_{\rm dep}\equiv{M}_{\rm g}/\langle\dot{M}_{\star}\rangle_{\Delta t}$. This estimate of $\tau_{\rm dep}$ can be related to the dynamics of individual gas parcels within the patch by noting that due to mass conservation, the instantaneous total mass of the gas parcels in the star-forming state, ${M}_{\rm sf}$, evolves as $\dot{M}_{\rm sf}=F_{\rm sf}-\dot{M}_{\star}$, where $F_{\rm sf}$ is the net instantaneous flux of the gas parcels through the star formation threshold and $\dot{M}_{\star}$ is the instantaneous SFR of all parcels. After averaging this expression over the time interval $\Delta t$, we get

[TABLE]

For brevity, in the following derivation, we omit explicit averaging, $\langle...\rangle_{\Delta t}$, but assume all fluxes and rates to be averaged over $\Delta t$.

In general, the average net flux of gas through the star formation threshold can be decomposed into positive and negative contributions,

[TABLE]

which correspond to the supply and removal of star-forming gas illustrated in Figure 1.

The positive flux $F_{\rm+}$ is controlled by a combination of global dynamical processes, e.g., gravitational instabilities, turbulence, cooling, etc., with a significant contribution from stellar feedback. The latter comes in the form of turbulence stirred by interactions of supernova-driven bubbles and by fountain outflows of gas from star-forming regions. The negative flux $F_{\rm-}=F_{\rm-,fb}+F_{\rm-,d}$ results from the destruction of star-forming regions both by feedback from young stars formed inside the regions, $F_{\rm-,fb}$, and by shearing due to large-scale turbulence or differential rotation, expansion of gas behind galactic spiral arms, and other dynamical processes not directly related to star formation inside the regions, $F_{\rm-,d}$. All fluxes can be parameterized with the characteristic timescales, i.e., $\tau_{+}$, $\tau_{\rm-,fb}$ and $\tau_{\rm-,d}$, on which gas is supplied to and removed from the star-forming state by the corresponding processes,

[TABLE]

where $f_{\rm sf}\equiv{M}_{\rm sf}/{M}_{\rm g}$ is the star-forming mass fraction.

To make the relation between star formation and stellar feedback explicit, we can also parameterize $F_{\rm-,fb}$ in a way similar to the parameterization of the mass outflow rate of feedback-driven galactic winds,

[TABLE]

where $\xi$ is the mass-loading factor and we used the definition of the average depletion time of star-forming gas, $\tau_{\star}\equiv{M}_{\rm sf}/\dot{M}_{\star}$. In the context of Equation (6) the mass-loading factor can also be interpreted as a relative rate of gas removal by feedback compared to the rate of star formation, i.e., $\xi\equiv\tau_{\star}/\tau_{\rm-,fb}$.

An imbalance between the net gas flux into the star-forming state, $F_{\rm sf}$, and the average SFR may result in the evolution of the star-forming mass, which we also parameterize with the characteristic timescale, $\tau_{\rm e,sf}$:

[TABLE]

The final expression for the global depletion time can be readily derived by substituting Equations (4-8) and $f_{\rm sf}=\tau_{\star}/\tau_{\rm dep}$ into Equation (3),

[TABLE]

where the sign in front of $\tau_{\star}/\tau_{\rm e,sf}$ reflects the sign of $\dot{M}_{\rm sf}$.

If we compare the terms in this equation with those in Equation (2) for the depletion time of a single gas parcel, $t_{\rm dep}=N_{\rm c}t_{\rm nsf}+\tau_{\star}$, their physical meaning becomes clear. The timescale $\tau_{+}$ is analogous to the time $t_{\rm nsf}$ that a gas parcel spends in the non-star-forming state, while the expression in parentheses is analogous to $N_{\rm c}=\tau_{\star}/t_{\rm sf}$, i.e., the average number of evolution cycles it would take for a single parcel to deplete its gas. Indeed, Equation (5) gives ${M}_{\rm nsf}=F_{\rm+}\tau_{+}$, which means that $\tau_{+}$ is the time over which all of the non-star-forming gas will reach the star-forming state. Thus, $\tau_{+}$ is analogous to the average $t_{\rm nsf}$ timescale for a collection of parcels. Likewise, the average rate at which the gas mass in the star-forming state is decreasing due to star formation, dispersal, and the overall evolution of the gas PDF during $\Delta t$ is given by ${M}_{\rm sf}/\tau_{\star}+{M}_{\rm sf}/\tau_{-}\pm{M}_{\rm sf}/\tau_{\rm e,sf}$, and the associated timescale $(1/\tau_{\star}+1/\tau_{-}\pm 1/\tau_{\rm e,sf})^{-1}$ corresponds to the average time that gas spends in this state. Thus, on average, gas will have to reach the star-forming state

[TABLE]

times, where we used Equation (6) and the definition of the mass-loading factor, $\xi\equiv\tau_{\star}/\tau_{\rm-,fb}$.

Equation (9) is the key expression of our framework. It states that the global depletion time is the sum of the total time that gas spends in the non-star-forming state over $N_{\rm c}$ cycles and the total time over which star-forming regions convert this gas into stars, $\tau_{\star}$.

This equation elucidates how long $\tau_{\rm dep}$ values can be reconciled with the relatively short local depletion times, $\tau_{\star}$, and even shorter dynamical timescales, $\tau_{+}$, discussed in the Introduction. The global depletion time is longer than the depletion time in star-forming regions, $\tau_{\star}$, due to the significant fraction of time that gas spends in the non-star-forming state. The global depletion time is longer than the timescale associated with dynamical processes supplying star-forming gas, $\tau_{+}$, because gas must evolve through the non-star-forming state $N_{\rm c}$ times, and $N_{\rm c}$ is large due to either efficient feedback, i.e., large $\xi$, or fast dynamical processes destroying star-forming regions, i.e., short $\tau_{\rm-,d}$ (see Equation 10).

When feedback dominates the removal of gas from the star-forming state, the number of cycles becomes $N_{\rm c}\sim\tau_{\star}/\tau_{\rm-,fb}$. This clarifies how feedback can self-regulate star formation, i.e., how $\tau_{\rm dep}$ can become insensitive to $\tau_{\star}$. Indeed, the timescale $\tau_{\rm-,fb}$ is proportional to the rate of energy and momentum injection by feedback, which, in turn, is set by the local rate of star formation, i.e., $\tau_{\star}$. Hence, $\tau_{\rm-,fb}\propto\tau_{\star}$, which renders $N_{\rm c}$ insensitive to $\tau_{\star}$. Thus, when $N_{\rm c}\tau_{+}\gg\tau_{\star}$, the depletion time, $\tau_{\rm dep}\approx N_{\rm c}\tau_{+}$, will be insensitive to $\tau_{\star}$.

In a nonequilibrium state, in which $\dot{M}_{\rm sf}>0$ ($<0$) during $\Delta t$, the term $\pm\tau_{\star}/\tau_{\rm e,sf}$ in Equation (9) accounts for the correction of the average rates estimated using the star-forming gas fraction, $f_{\rm sf}$, defined for the instantaneous masses ${M}_{\rm sf}$ and ${M}_{\rm g}$. This correction appears because, when $\dot{M}_{\rm sf}>0$ ($<0$), the actual average fraction of $\Delta t$ that gas spends in the star-forming state is smaller (higher) than $f_{\rm sf}$ and therefore more (fewer) transition cycles are required for depletion.

In a steady state, on the other hand, the gas distribution is stationary and the star formation rate is in equilibrium with the gas fluxes into and out of the star-forming state: $\dot{M}_{\rm sf}=F_{\rm sf}-\dot{M}_{\star}\approx 0$. In this case, $\tau_{\rm e,sf}\rightarrow\infty$ and the term $\tau_{\star}/\tau_{\rm e,sf}$ can be neglected in Equation (9). In such a steady state, $\tau_{\rm dep}={M}_{\rm g}/F_{\rm sf}$ (see Equation 3), and depletion time is determined by the net rate of gas inflow into the star-forming state, $F_{\rm sf}$. When $F_{\rm sf}$ is small, the depletion time is long. Galaxies as a whole reach the steady state with $\dot{M}_{\rm sf}\approx 0$ on the shortest of the timescales that control the global depletion time in Equation (9). Thus, globally, such an assumption is justified. However, individual ISM patches may deviate from the steady state, and the $\tau_{\star}/\tau_{\rm e,sf}$ term will be one of the sources of the scatter in depletion times.

In the remainder of the paper, we illustrate the framework described above using the results of an isolated galaxy simulation. Although the simulation adopts specific choices for many parameters, including resolution and prescriptions for star formation and feedback, the overall features and implications of our model do not depend on these specific choices. Our framework generically allows one to relate the depletion time on a large scale, e.g., the scale of an entire galaxy, to the star formation and feedback model that operates on a smaller scale, e.g., resolution scale of a simulation, where the distribution of local depletion times, $t_{\star}$, is defined.

3. Simulation

To illustrate the framework outlined above and elucidate the physical processes that give rise to long global depletion times, we use a simulation of an isolated $\sim L_{*}$-sized galaxy. We carried out the simulation with the adaptive mesh refinement $N$-body and gas dynamics code ART (Kravtsov.1999; Kravtsov.etal.2002; Rudd.etal.2008; Gnedin.Kravtsov.2011) and followed the evolution of an isolated gaseous disk in a live potential of a dark matter halo, stellar bulge, and stellar disk that are modeled with collisionless particles.

We adopt the initial conditions that were used in the AGORA code comparison project (agora2) and also in the studies of Agertz.etal.2013 and Semenov.etal.2016. Specifically, the isolated disk is initialized inside a dark matter halo with $v_{\rm c,200}=150{\rm~{}km~{}s^{-1}}$ and an initial concentration of $c=10$. The initial disk of old stars has an exponential density profile with a radial scale length of $r_{\rm d}\approx 3.4{\rm~{}kpc}$ and a vertical scale height of $h_{\rm d}=0.1r_{\rm d}$ with a total mass of $M_{\rm\star,d}\approx 3.4\times 10^{10}{\rm~{}M_{\odot}}$. The stellar bulge has an initial mass of $M_{\rm\star,b}\approx 4.3\times 10^{9}{\rm~{}M_{\odot}}$ that is distributed with a Hernquist density profile with $a=0.1r_{\rm d}$ (Hernquist.1990). The initial exponential gaseous disk has the same $r_{\rm d}$ and $h_{\rm d}$ as the stellar disk; its total mass is ${M}_{\rm g}\approx 8.6\times 10^{9}{\rm~{}M_{\odot}}$, which corresponds to the disk gas fraction of $f_{\rm g}\equiv{M}_{\rm g}/(M_{\rm\star,d}+{M}_{\rm g})=20\%$.

Gas evolution is governed by modified hydrodynamical equations that include terms related to cooling and heating, dynamical effects of subgrid turbulence, gas consumption by star formation, and injection of mass, momentum, and energy by feedback from young stars.

Cooling in the optically thin limit is implemented following the model of Gnedin.Hollon.2012. We assume a fixed metallicity of $Z=Z_{\odot}$ and constant background heating by interstellar radiation in the Lyman-Werner band with the photodissociation rate of $10^{-10}\ \rm s^{-1}$ (Stecher.Williams.1967). To model temperatures in dense self-shielded gas, we assume that extinction is proportional to the local column density of atomic gas, which we approximate as $nL_{\rm J,40}$, where $n$ is the gas number density in a cell and $L_{\rm J,40}$ is the local Jeans length with an applied temperature ceiling of $40{\rm~{}K}$ (model “L1a” in Safranek-Shrader.etal.2017).

In our simulation, we adopt the unresolved turbulence model of Schmidt.etal.2014 that dynamically follows the kinetic energy of gas motions on subgrid scales, $K$, as a separate hydrodynamical field akin to thermal energy. Subgrid turbulence is generated by random velocities on the resolution scale; it exerts pressure and viscous forces on resolved gas motions and decays on the local cell-crossing timescale, as motivated by simulations of both subsonic and supersonic turbulence (e.g., MacLow.etal.1998). For details about the model implementation in ART, we refer the reader to Semenov.etal.2016.

We adaptively resolve cells where the total gas mass exceeds $\sim 8\,300{\rm~{}M_{\odot}}$ and reach a maximum resolution of $\Delta=40{\rm~{}pc}$. Such a $\Delta$ is sufficient to resolve ISM structure down to densities of $n\sim 100-1\,000{\rm~{}cm^{-3}}$, and therefore we do resolve the dynamical processes that are sometimes claimed to limit the star-forming gas supply from the general ISM with average density of $n\sim 1{\rm~{}cm^{-3}}$. At the highest resolution level, we do not apply an artificial pressure floor in cold gas. Thus, the densities of star-forming regions are limited only by the effects of stellar feedback and the effective pressure due to thermal and both subgrid and resolved turbulent motions.

As a star formation prescription, we adopt the local rate ${\dot{\rho}}_{\star}=\rho/t_{\star}$, where $\rho$ is the gas density in a cell, and the local depletion time, $t_{\star}$, is related to the free-fall time, $t_{\rm ff}\equiv\sqrt{3\pi/32G\rho}$, with the efficiency per free-fall time, $t_{\star}=t_{\rm ff}/\epsilon_{\rm ff}$. Numerical and analytical models of star formation in turbulent GMCs generally predict a strong increase of $\epsilon_{\rm ff}$ with an increasing relative importance of gravity that facilitates star formation, in comparison to turbulence that provides support against gravity (for a review, see Padoan.etal.2014). For example, Padoan.etal.2012 found that the star formation efficiency of a turbulent cloud increases exponentially with a decreasing virial parameter, $\epsilon_{\rm ff}\approx\exp(-\sqrt{\alpha_{\rm vir}/0.53})$, where $\alpha_{\rm vir}$ is defined for a box with a side $\Delta$ as for a uniform sphere of radius $R=\Delta/2$:

[TABLE]

In our simulations, to apply the Padoan.etal.2012 fit in thermally supported gas, in the definition of $\alpha_{\rm vir}$ we consider contributions of both the sound speed, $c_{\rm s}$, and the explicitly modeled subgrid turbulent velocities, $\sigma_{\rm t}\equiv\sqrt{2K/\rho}$:

[TABLE]

Also, even though we are able to model $\epsilon_{\rm ff}$ following the Padoan.etal.2012 formula (see, e.g., Semenov.etal.2016; Li.etal.2017), in the simulation used here, we approximate the continuous exponential dependence of $\epsilon_{\rm ff}$ on $\alpha_{\rm vir}$ assuming a constant $\epsilon_{\rm ff}=1\%$ for $\alpha_{\rm vir}<\alpha_{\rm vir,sf}=10$, and $\epsilon_{\rm ff}=0$ elsewhere. This simulation was carried out as part of a simulation suite in which the $\epsilon_{\rm ff}$ value was varied systematically to explore its effect on star formation in galaxies. Using a constant $\epsilon_{\rm ff}$ and a sharp $\alpha_{\rm vir,sf}$ threshold makes the interpretation of simulation results easier, and we will present the results of this simulation suite in a forthcoming study. We explicitly checked that the global depletion times and the Kennicutt-Schmidt relations are similar in runs where $\epsilon_{\rm ff}$ follows the Padoan.etal.2012 fit and where we approximate this fit with a threshold. A qualitatively similar star formation prescription but with a different choice of parameters was studied by Hopkins.etal.2013.

We stress that the scenario of gas depletion described in Section 2 remains valid for any choice of star formation prescription, although in Section 5.2 we argue that such a prescription should be chosen carefully, as it is important for the prediction of realistic ISM properties. We note that our adopted threshold value, $\alpha_{\rm vir,sf}=10$, is consistent with the observed distribution of $\alpha_{\rm vir}$ in molecular clouds with sizes comparable to our resolution of $\Delta=40{\rm~{}pc}$: such clouds tend to have $\alpha_{\rm vir}\lesssim 10$ (MivilleDeschenes.etal.2016). The adopted value of $\epsilon_{\rm ff}=1\%$ is consistent with the average values deduced for observed actively star-forming clouds.

In our feedback prescription, we inject momentum and energy from the type II supernovae (SNe) at a uniform rate between 3 and 43 Myr after the formation of a stellar particle. The total number of SNe exploded over this time interval is computed assuming the Chabrier.2003 IMF. The thermal energy and radial momentum injected by each supernova are calibrated against simulations of a SN remnant evolution in a nonuniform medium (Martizzi.etal.2015), taking into account ambient gas density and with an additional boost of radial momentum by a fiducial factor of 5. We adopt such a boosting factor to compensate for the numerical loss of injected momentum and to account for the actual physical uncertainties of the radial momentum estimates (e.g., Gentry.etal.2017). In addition to SNe type II feedback, during evolution, stellar particles return a fraction of their mass following the prescription of Leitner.Kravtsov.2011.

Figure 2 shows the spatial distribution of gas number density, temperature, and subgrid turbulent velocity in our simulated galaxy at $t=500{\rm~{}Myr}$. The figure highlights the multiphase, dynamic structure of the ISM. Comparison with a simulation with no feedback, which we will present in a companion paper, shows that the structure of the ISM is significantly affected by stellar feedback. Its effect is manifested not only in the ubiquitous regions of hot, turbulent gas but also in the overall flocculent nature of the spiral pattern.

In our analysis, we estimate the molecular gas density in each cell using the KMT model (KMT1; KMT2; McKee.Krumholz.2010): $\rho_{\rm H_{2}}=\max[0,(1-0.75s/(1+0.25s))\rho]$, where at solar metallicity $s\approx 1.8/\tau_{\rm c}$ and $\tau_{\rm c}=320(\rho\Delta/{\rm g\ cm^{-2}})$. This model predicts a rapid increase of the molecular fraction at densities $n>10{\rm~{}cm^{-3}}$ indicated by the black contour in Figure 2. The resulting total mass fraction of molecular gas in our simulation is $f_{\rm H_{2}}\equiv M_{\rm H_{2}}/{M}_{\rm g}\sim 20\%$.

Subgrid turbulent velocities, $\sigma_{\rm t}$, range from $\lesssim 3{\rm~{}km~{}s^{-1}}$ in the diffuse ISM between the spiral arms to $\sim 30-300{\rm~{}km~{}s^{-1}}$ in hot SNe bubbles. In this simulation, supernovae do not explicitly inject turbulent energy, and high $\sigma_{\rm t}$ in hot bubbles are generated by the subgrid turbulence model. In the cold dense gas, turbulent velocities are supersonic and also vary significantly, $\sigma_{\rm t}\sim 5-15{\rm~{}km~{}s^{-1}}$. Strong subgrid turbulence in cold gas results in high values of $\alpha_{\rm vir}$ and, according to our star formation criterion, $\alpha_{\rm vir}<\alpha_{\rm vir,sf}=10$, only $\sim 40\%$ of all molecular gas mass is star-forming at any given moment. Such star-forming regions are shown in the left panel of Figure 2 with green contours.

The total SFR of our model galaxy is $\dot{M}_{\star}\sim 1-2{\rm~{}M_{\odot}~{}yr^{-1}}$, which translates to global depletion times for the total and molecular gas of $\tau_{\rm dep}\sim 4-8{\rm~{}Gyr}$ and $\tau_{\rm dep,H_{2}}\sim 1-2{\rm~{}Gyr}$. The values of the depletion times and molecular fraction are in the ballpark of the typical values observed in nearby spiral galaxies (e.g., wong_blitz02; Bigiel.etal.2008; Leroy.etal.2013).

Figure 3 compares the Kennicutt-Schmidt relation between the surface densities of the SFR and H I + H2 and H2 gas in our simulation to the observed relations in the Milky Way and nearby spiral galaxies. Our results are in a good agreement in both normalization (i.e., the depletion time value) and slope. Note, in particular, that the linear relation between $\dot{\Sigma}_{\star}$ and $\Sigma_{\rm H_{2}}$ emerges from the nonlinear star formation prescription adopted in our simulation: ${\dot{\rho}}_{\star}\propto\rho^{1.5}$. In Section 4.4, we consider the origin of the linear relation and show that it results from the particular behavior of the gas density distribution shaped by stellar feedback.

For a consistent comparison with the observed KSR for H I + H2 gas, in our simulation, we defined neutral hydrogen to be all nonmolecular gas denser than $n_{\rm H,SSh}$, given by Equation (13) in Rahmati.etal.2013. This threshold corresponds to the gas self-shielded from the far ultraviolet (FUV) background with the adopted photoionization rate $\Gamma=10^{-10}\ \rm s^{-1}$. We also excluded all neutral hydrogen that is colder than $1000{\rm~{}K}$ assuming that it constitutes the optically thick cold neutral medium (CNM) not included into the observed measurements of $\Sigma_{\rm HI}$. Our temperature threshold is somewhat higher than the CNM temperature estimated in real galaxies ($\lesssim 300{\rm~{}K}$, e.g., Wolfire.etal.2003) because, in our simulation, we do not resolve the transition between warm and cold neutral gas phases, which results in intermediate gas temperatures on the resolution scale. The particular value of the temperature threshold was chosen to select $\sim 40\%$ of the neutral hydrogen mass, which is close to the CNM mass fractions estimated in the Milky Way and nearby galaxies (e.g., Heiles.Troland.2003; Braun.2012; Pineda.etal.2013; Sofue.2017).

In the analyses presented below, we consider the processes in the ISM during the time interval between $400$ and $600{\rm~{}Myr}$, which is short enough to neglect the effect of global gas consumption on the ISM dynamics and long enough to accumulate sufficient statistics of rare events and make the considered distributions representative and well-sampled.

To study the detailed dynamics of individual gas parcels, we use gas tracer particles that are passively advected with the local flow velocity interpolated to the positions of the particles with the cloud-in-cell scheme. We populate the disk with $10^{5}$ tracer particles uniformly initialized within $R<8{\rm~{}kpc}$ after $300{\rm~{}Myr}$ of disk evolution when the transients related to the initial off-equilibrium state had dissipated away. After initialization, we wait for $100{\rm~{}Myr}$ to let the tracers equilibrate with the gas density distribution. At that point, the distributions of tracer densities in radial annuli resemble the gas density PDF in computational cells, and the distribution of tracer properties thus can be considered to be a good approximation of the mass-weighted PDF of gas properties.

We average the distribution of tracers to construct statistics, such as their PDF and fluxes in the $n$ – $\sigma_{\rm tot}$ phase diagram, between $400$ and $600{\rm~{}Myr}$ with $1{\rm~{}Myr}$ steps. We checked that at every moment between $400$ and $600{\rm~{}Myr}$ phase distributions of gas and tracer particles resembled their averaged versions, which means that the galaxy remained in approximate equilibrium over the considered period of time.

We focus on the evolution of gas in the phase plane of density, $n$, and total velocity dispersion, $\sigma_{\rm tot}$, because the position of a gas parcel in this plane determines its internal consumption time, $t_{\star}$, according to our star formation prescription. To quantify gas motions in the $n$ – $\sigma_{\rm tot}$ plane, we measure the derivatives $d\log n/dt$ and $d\log\sigma_{\rm tot}/dt$ probed by each tracer every $1{\rm~{}Myr}$. To estimate the average local flow rates of gas, we accumulate fluxes corresponding to these derivatives and normalize them by the local density of tracers in the $n$ – $\sigma_{\rm tot}$ plane.

To characterize actual fluxes that supply and remove star-forming gas, in addition to the total flux of tracers, we separately track the fluxes of tracers with decreasing or increasing $\alpha_{\rm vir}$. We quantify the magnitudes of these fluxes with the characteristic evolution timescale, on which $\alpha_{\rm vir}$ changes by an order of magnitude at a given rate,

[TABLE]

where we average the derivative of $\log_{10}\alpha_{\rm vir}$ taking into account only tracers with decreasing (increasing) $\alpha_{\rm vir}$ to compute $\tau_{\alpha,<0}$ ($\tau_{\alpha,>0}$).

4. Results

4.1. The origin of long global depletion times

As we discussed at the end of Section 2, in steady state, when $\langle\dot{M}_{\rm sf}\rangle_{\Delta t}\approx 0$ on the timescale $\Delta t$ over which the SFR is estimated, long global depletion times, $\tau_{\rm dep}\equiv{M}_{\rm g}/\dot{M}_{\star}={M}_{\rm g}/F_{\rm sf}$, originate from a small net flux of gas into the star-forming state, $F_{\rm sf}$. In principle, $F_{\rm sf}$ could be small if the rate at which gas evolves toward the star-forming state were set by a slow “bottleneck” process. However, as Figure 4 shows, in simulations with efficient feedback, gas rapidly transitions between the star-forming and non-star-forming states, and a small $F_{\rm sf}$ results from a near-cancellation of large opposite fluxes into and out of the star-forming state.

In this figure, we plot the distribution of gas tracer particles within the disk in the plane of gas number density, $n$, and $\sigma_{\rm tot}=\sqrt{\sigma_{\rm t}^{2}+c_{\rm s}^{2}}$, that can be viewed as an effective temperature including both thermal and turbulent gas motions on subgrid scales. The gas distribution spans a wide range of densities, $\sigma_{\rm tot}$, and temperatures and has two distinct peaks. The peak at low densities, $n\sim 1{\rm~{}cm^{-3}}$, corresponds to diffuse, warm, subsonic ($\sigma_{\rm t}\lesssim c_{\rm s}$) gas at temperature $T\sim 10^{4}{\rm~{}K}$. The gas in the second peak at $n>10{\rm~{}cm^{-3}}$, on the other hand, is cold ($T\lesssim 100$ K) and supersonic ($\sigma_{\rm t}>c_{\rm s}$).

According to our star formation prescription, the star-forming gas has $\alpha_{\rm vir}<\alpha_{\rm vir,sf}=10$. Such gas in Figure 4 resides below the thick green dotted line. The net mass flux of gas in the $n-\sigma_{\rm tot}$ plane is visualized by the arrows in the top panel of Figure 4, where the length of the arrows is equal to the distance tracers would traverse in $5{\rm~{}Myr}$ for a given flux. The figure shows that arrows are rather small throughout most of the phase space occupied by tracers and are particularly small near the thick green dotted line. This means that the net evolution of gas in the $n$ – $\sigma_{\rm tot}$ plane is slow and the net flux through the star formation threshold, $F_{\rm sf}$, is small. This small net flux results in the long global depletion timescales exhibited by our simulated galaxy, $\tau_{\rm dep}\sim 5{\rm~{}Gyr}$ and $\tau_{\rm dep,H_{2}}\sim 1{\rm~{}Gyr}$ (see Section 3 and Figure 3).

However, the middle and bottom panels of Figure 4 show that the small net $F_{\rm sf}$ results from the near-cancellation of two opposite fluxes. These panels show the fluxes of only those tracers in which $\alpha_{\rm vir}$ is decreasing, $F_{\rm+}$, or increasing, $F_{\rm-}$, and these fluxes are significantly stronger than the net flux in the top panel. A typical tracer evolves toward and away from the star-forming state on a timescale of order $\tau_{\alpha,\gtrless 0}\sim 5-30$ Myr, consistent with the estimates of the timescales of relevant processes in Section 1. Thus, the rate of gas supply from the diffuse, warm ISM to the star-forming state cannot be the factor limiting the global star formation rate, as envisioned by Saitoh.etal.2008. Instead, gas generally evolves from the diffuse to the star-forming state on a timescale of tens of Myr, much shorter than the global depletion time. The latter is long because gas rapidly leaves the star-forming state at the rate that nearly cancels the rate at which gas is reaching this state. In the next section, we consider the processes that drive the fast gas evolution in more detail.

4.2. Dynamical processes shaping ISM

The average gas flow patterns shown with arrows in Figure 4 result from the statistical averaging of the complicated trajectories of individual tracer particles. The particular shapes of such trajectories vary depending on local conditions and specific physical processes that govern gas evolution.

In our simulation, gas evolution between diffuse, warm, subsonic and dense, cold, supersonic ISM phases is governed by large-scale disk instabilities and the turbulent flows generated by them. The evolution of gas in the dense phase is determined predominantly by stellar feedback that disperses star-forming regions, drives large-scale ISM turbulence, and launches fountain-like outflows.

In the following subsections, we consider these processes using three illustrative tracer trajectories integrated over $200{\rm~{}Myr}$ and shown in Figure 5. We chose these particular tracers because their evolution over the considered period of time is governed predominantly by the same process over several consequential cycles of compression and expansion.

4.2.1 Compression and expansion due to spiral arms

The blue line in the top left panel of Figure 5 shows an example of the $\alpha_{\rm vir}$ evolution followed by a tracer that swings between the subsonic and supersonic phases during cycles of compression and expansion as it enters and exits the spiral arms. The trajectory of this tracer in the $n$ – $\sigma_{\rm tot}$ plane is shown with the same color in the right panel.

Due to strong compression, gas entering a spiral arm rapidly cools down and loses the thermal support that initially dominates in the subsonic ISM. At the same time, initially low subgrid turbulent velocities of the subsonic ISM, $\sigma_{\rm t}\lesssim 3{\rm~{}km~{}s^{-1}}$, rapidly grow due to compressional heating (see the detailed discussion in Section 3.1 of Semenov.etal.2016 and also Robertson.Goldreich.2012). At $n\sim 10{\rm~{}cm^{-3}}$, when subgrid turbulent velocities become comparable to the thermal speed, gas detaches from the lower envelope of the distribution shown in Figure 4 and enters the supersonic ISM phase. Similarly, when gas leaves a spiral arm, it expands, and subgrid turbulent velocities decrease. Eventually, under the influence of expansion and interstellar FUV heating, gas returns to the subsonic ISM phase with $n\sim 1{\rm~{}cm^{-3}}$ and $T\sim 10^{4}{\rm~{}K}$.

The actual transition of gas between the subsonic and supersonic phases is fast, as it is controlled by a strong compression rate in the spiral arms and short cooling times at $n>1{\rm~{}cm^{-3}}$. Hence, the rate at which diffuse gas is promoted into the dense phase is mostly determined by the time that gas waits between subsequent passages of the spiral arms,

[TABLE]

where $V_{\rm gas}\equiv v_{\rm gas}-v_{\rm pat}$ is the speed of gas relative to the spiral waves pattern, and we set $m=6$, as our simulated galaxy develops six spiral arms.

The typical time that gas spends inside a spiral arm before expansion contributes to the dynamical rate of gas removal from the star-forming state, $\tau_{\rm-,d}$. This timescale depends on the spiral arm width, gas velocity, and the angle at which gas flows inside the arm. Depending on local conditions, this timescale can be as long as a few tens of Myr.

In Figure 6 we plot distributions of $n$ and $\sigma_{\rm tot}$ separately for tracers residing at different galactic radii and therefore experiencing different ISM conditions. The distribution in the outer disk (top panel) is shaped predominantly by the compression and expansion due to the spiral arms. Specifically, most of the gas mass in the outer disk resides in the diffuse subsonic phase and forms a peak at $n\sim 1{\rm~{}cm^{-3}}$ and $T\sim 10^{4}{\rm~{}K}$. The tail extending along the lower envelope of the distribution toward the dense supersonic phase corresponds to the gas currently being compressed in the spiral arms. As the figure also shows, the compression of diffuse gas in the spiral arms is only relevant at large radii, whereas closer to the disk center, less gas remains in the diffuse phase, and this process becomes much less important.

4.2.2 SNe-induced shocks and ISM turbulence

We find that the evolution of dense, supersonic gas in the $n$ – $\sigma_{\rm tot}$ plane is dominated by the turbulence that is driven by stellar feedback. Injection of momentum by SNe in a star-forming region results in a rapid expansion of gas until the region is eventually dispersed. Shocks associated with expanding bubbles compress gas in the disk plane, which may induce new episodes of star formation and subsequent SN explosions. The turbulence resulting from overlapping and interacting bubbles makes gas parcels oscillate in fast cycles, as illustrated by the green trajectory in Figure 5.

The characteristic timescale between subsequent compressions of ISM gas by such expanding SN shocks corresponds to

[TABLE]

where $L$ is a typical separation between bubbles (see, e.g., the temperature map in the middle panel of Figure 2) and $v_{\rm shell}$ is a typical velocity of shells on a scale $L$.

Compression and expansion of gas in the turbulent ISM is accompanied by the increase and decrease of turbulent velocity dispersion. As a result, averaging of such large-scale turbulent motions over many tracers results in a prominent clockwise whirl-like pattern of arrows around the peak of the PDF in the cold (blue) part of the diagram (see Figures 4 and 6). Closer to the peak center, the net flux magnitude decreases due to the averaging between fast motions of many tracers at the different stages of their turbulent compression-expansion cycles.

As Figure 6 shows, such a whirlwind pattern is most prominent at $R\leq 6$ kpc. Thus, the feedback-driven turbulence and associated compression and expansion of gas are dominant processes at these radii in the cold, supersonic gas. The ISM at these radii has a complex structure (see Figure 2) reflecting the chaotic turbulent nature of the gas.

4.2.3 Feedback-driven fountain outflows

Supernova feedback also affects some of the gas by accelerating it in the direction perpendicular to the disk plane. Such gas expands in fountain-like outflows but eventually cycles back to the ISM under the influence of the disk potential. Interactions of such outflows with the halo gas adjacent to the disk result in an increase of small-scale turbulent velocities that quickly dissipate when the gas falls back onto the disk.

An example of a tracer trajectory during expansion and subsequent recycling of a fountain outflow is shown with the red lines in Figure 5. This particular tracer was ejected and recycled twice, at $\sim 20-80$ and $\sim 130-150{\rm~{}Myr}$. In each event, after its star-forming region was dispersed by feedback, this tracer acquired a moderate vertical velocity of $v_{z}\sim 50{\rm~{}km~{}s^{-1}}$ and elevated as high as $\sim 400{\rm~{}pc}$ above the disk plane, i.e., a few scale heights, before falling back onto the disk. At the highest elevation point, the gas in these outflows expands only to the densities comparable to those of the diffuse subsonic ISM phase, $n\sim 0.2-2{\rm~{}cm^{-3}}$ and its virial parameter reaches the values of $\alpha_{\rm vir}\sim 10^{3}-10^{4}$ due to the strong turbulence generated by the interaction of the expanding outflow with the surrounding gas. Outflows launched by feedback from regions of more vigorous star formation reach even lower $n$ and higher $\alpha_{\rm vir}$.

The timescale of the fountain cycle can be estimated as a dynamical time in the gravitational field of a massive infinite sheet of constant surface density $\Sigma_{\rm tot}$, corresponding to the local total surface density of the disk,

[TABLE]

where $v_{z}$ is the initial vertical velocity of gas in the outflow and $\Sigma_{\rm tot}=\Sigma_{\star}+\Sigma_{\rm g}\sim 200{\rm~{}M_{\odot}}{\rm~{}pc}^{-2}$ is the typical total surface density of gas and stars in our simulated galaxies.

The averaging of trajectories between many gas parcels constituting fountain-like outflows results in a tail of the distribution directed from the star-forming state toward the lower densities and higher $\sigma_{\rm tot}$. The total flux of tracers forms a prominent counterclockwise vortex inside this tail that is clearly seen in Figures 4 and 6.

Figure 6 shows that at all radii within the disk, some fraction of gas evolves in the manner discussed above, which indicates the existence of fountain-like outflows. At larger radii, where the SFR is slower, the outflows are less prominent but still visible as a net flux of tracers directed toward lower densities along the top envelope of the distribution shown in the upper panel. Gas in such outflows at large radii usually returns to the diffuse, warm, subsonic ISM between the spiral arms. Closer to the center, outflows are ubiquitous, and, after falling back, their gas directly rejoins the tumultuous large-scale turbulent motions of dense, supersonic gas.

4.3. Duration and number of evolution cycles

Typical tracer trajectories considered in the previous section explicitly confirm that during the evolution, gas parcels perform many fast cycles and rapidly explore a significant portion of the PDF extent, frequently switching between non-star-forming and star-forming states. As we discussed in Section 2, the distribution of the time that the gas parcels spend on each cycle in these states, $t_{\rm nsf}$ and $t_{\rm sf}$, determines the global depletion time of the galaxy.

In Figure 7 we plot the distribution of $t_{\rm nsf}$ and $t_{\rm sf}$ directly measured from the trajectories of all tracers as the time between consequential crossings of the star formation threshold. The results in the previous section indicate that the mix of the processes governing gas evolution may change with the galactic radius, $R$, and surface density, $\Sigma_{\rm g}$. Thus, to explore possible trends, we plot the distributions of timescales as a function of $R$ and $\Sigma_{\rm g}$.

The distribution of $t_{\rm nsf}$ shown in the top panel indicates that the majority of tracers spend less than $100{\rm~{}Myr}$ in the non-star-forming stage of evolution during each cycle. At higher average surface densities closer to the disk center, this time is even shorter, $t_{\rm nsf}\lesssim 50{\rm~{}Myr}$, with a very low median value (thick red line).

Colors in the top panel show the average fraction of time that gas tracers spend in the dense phase, $n>10{\rm~{}cm^{-3}}$, over the non-star-forming stage of evolution. A blue color at small radii implies that gas preferentially stays in the dense, molecular phase even when it does not form stars. This is also evident from the bottom panel of Figure 6, which shows that only a small fraction of gas expands to $n<10{\rm~{}cm^{-3}}$ and it does so as a part of fountain outflows.

At larger radii, the relatively slow rate of star-forming gas replenishment via compression in the spiral arms becomes important, and the median $t_{\rm nsf}$ increases to $\sim 80{\rm~{}Myr}$. Gas governed by this process spends significant time in the diffuse subsonic ISM, and such tracer trajectories occupy the areas of the longest $t_{\rm nsf}$ at $R>4$ kpc (green color in Figure 7). However, as indicated by the blue color, at such radii many tracers still perform short cycles with $t_{\rm nsf}<50{\rm~{}Myr}$ without leaving the dense phase.

The increase of the $t_{\rm nsf}$ median value is consistent with the scaling proportional to the free-fall time at the mean or midplane density at a given radius, $t_{\rm ff,0}\propto\rho_{0}^{-1/2}=(\Sigma_{\rm g}/2h_{\rm d})^{-1/2}$, shown by the dotted line in Figure 7 (see also Saitoh.etal.2008). Such scaling is sometimes adopted in analytical models of galactic star formation to define the timescale on which star-forming regions are created (e.g., Krumholz.etal.2012; Elmegreen.2015). As we discussed above, $t_{\rm nsf}$ in our simulations is set by both stellar feedback that drives turbulence and dynamical processes within the ISM. The scaling of the median $t_{\rm nsf}$ with density indicates that gravity and the associated timescale plays at least some role in setting the time that gas spends in the non-star-forming state. For example, the fall of the gas driven out in a fountain outflow back to the disk will occur on a timescale of order $\sim t_{\rm ff,0}$.

The bottom panel of Figure 7 shows the distribution of time spent by tracers in the star-forming state on each cycle. This timescale is close to the typical “lifetime” of star-forming regions and is quite short: $t_{\rm sf}\lesssim 20{\rm~{}Myr}$ or $2-4$ free-fall times at the typical densities of star-forming regions. The fact that $t_{\rm sf}$ is, on average, significantly shorter than $t_{\rm nsf}$ is consistent with the small mass fraction of star-forming gas.

As we discussed in Section 2, the average time that a gas parcel spends in the star-forming state on a single cycle determines the total number of such cycles required for complete depletion as $N_{\rm c}=\tau_{\star}/t_{\rm sf}$. For every tracer on each passage through the star-forming stage, we estimate this number as the inverse fraction of mass depleted during the passage,

[TABLE]

where the integral is accumulated for each tracer particle while it is in the star-forming state between subsequent crossings of the star formation threshold. The resulting distribution of $N_{\rm c}$ is shown by the colors in the bottom panel of the figure.

In agreement with our model, typical $N_{\rm c}\sim 50$ and the lifetimes of gas in the star-forming state, $t_{\rm sf}\sim 10-20{\rm~{}Myr}$, are consistent with the range of the star-forming gas depletion times, $\tau_{\star}\sim N_{\rm c}t_{\rm sf}\sim 300-500$ Myr, obtained in our simulation (see the next subsection). In addition, assuming $t_{\rm nsf}\sim 50-100{\rm~{}Myr}$, Equation (2) for the typical depletion time of a gas parcel gives a value of $\sim 2-5{\rm~{}Gyr}$, which is consistent with the actual global depletion time obtained in our simulation (see Section 3).

4.4. Emergence of linear $\dot{\Sigma}_{\star}-\Sigma_{\rm H_{2}}$ relation

One of the most intriguing results of our simulation is the emergence of the linear $\dot{\Sigma}_{\star}-\Sigma_{\rm H_{2}}$ relation consistent with observations (see the right panel of Figure 3), even though on small scales the star formation rate is assumed to scale nonlinearly with the gas density: ${\dot{\rho}}_{\star}=\epsilon_{\rm ff}\rho/t_{\rm ff}\propto\rho^{1.5}$. This result provides a counterexample to the arguments that in simulations, the slope of the Kennicutt-Schmidt relation on large and small scales should be the same (Schaye.DallaVecchia.2008; Gnedin.etal.2014).

The observed linearity of the molecular KSR is often explained by the “counting argument,” which posits that all molecular gas is in star-forming clouds of approximately the same depletion time. In this case, the surface density of the SFR will vary linearly with the surface density of molecular gas, as long as the geometric covering fraction of star-forming clouds is less than unity.

However, if not all of the molecular gas is star-forming, as is the case in our simulations and is likely the case in observed galaxies, the explanation for the linearity must be more nuanced. Analogously to Equation (1), the depletion time of molecular hydrogen is given by

[TABLE]

where $f_{\rm sf,H_{2}}\equiv\Sigma_{\rm sf}/\Sigma_{\rm H_{2}}$. The linearity of the molecular KSR, or, equivalently, the independence of the $\tau_{\rm dep,H_{2}}$ of $\Sigma_{\rm H_{2}}$, must therefore arise from a lack of dependence of $\tau_{\star}$ and $f_{\rm sf,H_{2}}$ on $\Sigma_{\rm H_{2}}$ or a cancellation of any such dependence in their ratio.

The depletion time of star-forming gas is given by the average over the inverse local depletion time distribution in star-forming regions, $\tau_{\star}\equiv\Sigma_{\rm sf}/\dot{\Sigma}_{\star}=\langle 1/t_{\star}\rangle^{-1}_{\rm sf}$. Thus, if the distribution of $t_{\star}$ is independent of $\Sigma_{\rm H_{2}}$, so is $\tau_{\star}$, regardless of the actual shape of the distribution. In our simulation, $t_{\star}=t_{\rm ff}/\epsilon_{\rm ff}$, where $\epsilon_{\rm ff}=\mathrm{const}$, and thus the distribution of $t_{\star}\propto n^{-1/2}$ is set by the density distribution of star-forming gas.

The PDFs of the star-forming gas density in radial annuli with widths of 500 pc are shown in Figure 8. The figure shows that $\Sigma_{\rm H_{2}}$ spans two orders of magnitude, while the average star-forming gas density and the shape of the PDF do vary, but the variation is quite mild. Thus, $\tau_{\star}$ depends on $\Sigma_{\rm H_{2}}$, but weakly, as shown by the thin line in Figure 9.

This dependence is weak for two reasons. First, the average density of star-forming gas in all of the ISM environments in our simulation is limited by the stellar feedback (at high densities) and FUV background (at $n<20{\rm~{}cm^{-3}}$) to a nearly constant range of $n\sim 10-10^{3}\ \rm cm^{-3}$, comparable to the typical densities of observed star-forming regions. Second, as in the counting argument, $\Sigma_{\rm H_{2}}$ varies mostly due to the variation of the geometric covering fraction of molecular gas, not due to the change of its density.

As the dotted line in Figure 9 shows, the star-forming mass fraction, $f_{\rm sf,H_{2}}$, also varies weakly with $\Sigma_{\rm H_{2}}$, with a dependence on $\Sigma_{\rm H_{2}}$ similar to that of $\tau_{\star}$. In Equation (18) the weak trends of $\tau_{\star}$ and $f_{\rm sf,H_{2}}$ with $\Sigma_{\rm H_{2}}$ cancel, resulting in a constant $\tau_{\rm dep,H_{2}}$, which explains the origin of the linear $\dot{\Sigma}_{\star}-\Sigma_{\rm H_{2}}$ relation.

This almost exact cancellation of the $f_{\rm sf,H_{2}}$ and $\tau_{\star}$ trends implies that $f_{\rm sf,H_{2}}\propto\tau_{\star}$. Some covariance of $f_{\rm sf,H_{2}}$ with $\tau_{\star}$ is expected because $\tau_{\star}$ controls the rate of feedback energy and momentum injection. As a result, longer $\tau_{\star}$ results in less turbulence, i.e., smaller $\alpha_{\rm vir}$ and larger $f_{\rm sf,H_{2}}$. This qualitative explanation, however, will need to be investigated and tested more thoroughly with simulations and will be a subject of our future study.

5. Discussion

The short, $\sim 10-100$ Myr, timescales of the physical processes driving the evolution of gas in the ISM (see the top panel of Figure 7) indicate that the ISM is vigorously “boiling” when considered on the global depletion timescale. During this vigorous evolution, gas cycles between non-star-forming and star-forming stages and spends only $t_{\rm sf}\sim 5-15{\rm~{}Myr}$ in the star-forming stage on each cycle (see the bottom panel of Figure 7), which is consistent with the short lifetimes derived for observed GMCs (e.g., Kawamura.etal.2009; Murray.2011; Schruba.etal.2017).

Observational estimates of the integral star-formation efficiency during a star-forming stage, defined for a given star-forming region containing a gas mass of $m_{\rm g}$ and a mass of formed young stars of $m_{\star}$ as $\epsilon_{\rm int}\equiv m_{\star}/(m_{\rm g}+m_{\star})$, give $\epsilon_{\rm int}\sim 1-20\%$ (e.g., Evans.etal.2009; Lada.etal.2010). This fraction is even smaller in less-efficient clouds (e.g., Rebolledo.etal.2015; Lee.etal.2016; Vutisalchavakul.etal.2016). A similar range of $\epsilon_{\rm int}$ is also obtained in simulations of star cluster formation (Gavagnin.etal.2017), models of star formation in GMCs (e.g., Zamora.Vazquez.Semadeni.2014), and cosmological simulations of a Milky Way-sized galaxy that resolve the growth of globular clusters and self-consistently capture its termination by stellar feedback (Li.etal.2017).

Such values of $\epsilon_{\rm int}$ imply that gas parcels must undergo $N_{\rm c}\sim\epsilon_{\rm int}^{-1}\sim 5-100$ cycles transitioning from the non-star-forming to star-forming state before they convert their gas into stars. This number of cycles is also consistent with the typical depletion times of star-forming gas, $\tau_{\star}\equiv\langle 1/t_{\star}\rangle^{-1}_{\rm sf}$, and lifetimes, $t_{\rm sf}$, derived for observed star-forming regions, $N_{\rm c}=\tau_{\star}/t_{\rm sf}$. Specifically, depletion times of gas in observed star-forming regions are estimated to be $t_{\star}\sim 50-500{\rm~{}Myr}$ (e.g., Evans.etal.2009; Evans.etal.2014; Lada.etal.2010; Lada.etal.2012; Heiderman.etal.2010; Gutermuth.etal.2011; Schruba.etal.2017); for $t_{\rm sf}\sim 5-15{\rm~{}Myr}$ quoted above, these $t_{\star}$ give $N_{\rm c}\sim(50-500)/(5-15)\sim 3-100$. Such $N_{\rm c}$ are in the ballpark of the $N_{\rm c}$ range that we estimate for the gas in our simulations (see the bottom panel of Figure 7), although we note that the number of cycles in our simulation can be somewhat overestimated due to the artificially sharp threshold in the definition of star-forming gas.

Note that specific values of $t_{\rm nsf}$, $t_{\rm sf}$ and $N_{\rm c}$ depend on the scale, $l$, on which the small-scale distribution of $t_{\star}$ is defined. Clearly, if we consider the evolution of gas parcels on the scale of protostellar cores, $\sim 0.01{\rm~{}pc}$, the star-forming stage of evolution will correspond to the formation of one or a handful of stars, which will consume most of the gas in a single event. The gas parcels on this scale will spend a long time in the non-star-forming stage and will consume their gas in one or a few cycles, $N_{\rm c}=1+\xi$, where $\xi$ corresponds to the ratio of the protostellar core mass returned back to the ISM to the mass of the formed star.

The key expression of our model (Equation 9) was derived by applying the mass conservation law to all star-forming gas in a galaxy or in a given ISM patch. Mass conservation can also be invoked to build a model for the formation, evolution, and destruction of individual GMCs (e.g., Feldmann.Gnedin.2011; Zamora.etal.2012; Lee.etal.2016).

Our model is based on mass conservation alone, and thus the overall physical explanation for long gas depletion times presented in Sections 2 and 4 does not require the assumption of dynamical equilibrium. Indeed, when a gas parcel undergoes some inherently nonequilibrium process, such as compression in a spiral arm, the parcel’s depletion time will still be given by $t_{\rm dep}=(\tau_{\star}/t_{\rm sf})(t_{\rm nsf}+t_{\rm sf})=\tau_{\star}(t_{\rm nsf}/t_{\rm sf}+1)$, and therefore $t_{\rm dep}$ will be long if $t_{\rm nsf}\gg t_{\rm sf}$ and/or $\tau_{\star}$ is long. The absence of a long-term equilibrium requirement is an essential difference of our framework from the models for the Kennicutt-Schmidt relation that rely on the assumption of self-regulation to the vertical or Toomre.1964 equilibrium state (e.g., Ostriker.Shetty.2011; FaucherGiguere.etal.2013).

In fact, our model explicitly accounts for the deviations from the equilibrium state in which $\dot{M}_{\rm sf}\approx 0$. Such deviations, along with fluctuations of other quantities that enter Equation (9), can be important sources of the depletion time scatter. These deviations can be substantial for individual ISM patches, which generally will not be in equilibrium, even if a galaxy as a whole is. Additional scatter can arise due to observational tracers sampling different stages of gas evolution incompletely (Kruijssen.Longmore.2014).

We note, however, that unlike the models of Ostriker.Shetty.2011 and FaucherGiguere.etal.2013, our conceptual framework cannot quantitatively predict the depletion time by itself. It only elucidates how the depletion time is related to the timescales of the processes driving gas evolution. The variables through which this relation is parameterized can be either calibrated in simulations, as is done in this study, or be derived in analytical models (e.g., Zamora.etal.2012). Nevertheless, as we discuss in the following subsections, our framework is very useful for interpreting and explaining a number of puzzling facts about star formation in both observed and simulated galaxies.

5.1. Implications for observations

Rapid cycling of gas between non-star-forming and star-forming states explains the large discrepancy between long global depletion times of $\gtrsim 1{\rm~{}Gyr}$ and short, $\sim 10-100{\rm~{}Myr}$, timescales associated with the dynamical processes in the ISM. Only a small fraction of gas is converted into stars during each cycle, and therefore gas would have to go through a large number of cycles to be depleted.

Our model also naturally explains the difference between observed local depletion times of (mostly molecular) gas in star-forming regions, $t_{\star}\sim 50-500{\rm~{}Myr}$, and global depletion times of both total gas, $\tau_{\rm dep}\sim 2-10{\rm~{}Gyr}$, and molecular gas, $\tau_{\rm dep,H_{2}}\sim 1-3{\rm~{}Gyr}$. The global depletion times, $\tau_{\rm dep}$ and $\tau_{\rm dep,H_{2}}$, are longer than the average gas depletion time in star-forming regions, $\tau_{\star}$, due to the significant fraction of time that gas spends in the non-star-forming state (see Equation 2). This implies that only a fraction of total and molecular gas is forming stars at any given moment. For example, the observed values of $t_{\star}$ and $\tau_{\rm dep,H_{2}}$ indicate that only $f_{\rm sf,H_{2}}\equiv{M}_{\rm sf}/M_{\rm H_{2}}=\tau_{\star}/\tau_{\rm dep,H_{2}}\sim 5\%-50\%$ of molecular gas is forming stars. The range of $f_{\rm sf,H_{2}}$ in our simulation is consistent with this estimate (see the dotted line in Figure 9), with the non-star-forming molecular gas state corresponding to strongly turbulent cold gas. Diffuse, non-star-forming molecular gas is also observed in the Milky Way, but its mass fraction ($\sim 25\%$ from RomanDuval.etal.2016) is a factor of $\sim 2-4$ too low to explain the discrepancy between global and local depletion times of molecular gas. Thus, a substantial fraction of non-star-forming molecular gas must be missed in such observations.

As pointed out by Kruijssen.Longmore.2014, a model considering different evolutionary stages and corresponding chemical phases of the ISM gas can be used to interpret the dependence of the depletion time and its scatter on the averaging scale (see, e.g., Schruba.etal.2017). This dependence in observed galaxies can also be used as a stringent test of the star formation and feedback implementation in galaxy formation models.

Our model for gas depletion time provides a natural framework for predicting and interpreting trends with galaxy properties and redshift. For instance, we show that the duration of cycles decreases with increasing surface density as $t_{\rm nsf}\propto\Sigma_{\rm g}^{-0.5}$, which is accompanied by a milder but nevertheless nonnegligible decrease in $\tau_{\star}$ (see Section 4.4 and the thin line in Figure 9). This means that the observed decrease in global depletion times in high-redshift and starburst galaxies (e.g., Kennicutt.1998; Bouche.etal.2007; Genzel.etal.2010; Tacconi.etal.2017) can be explained by shorter dynamical timescales, $t_{\rm nsf}$, and star-forming gas depletion times, $\tau_{\star}$, associated with high-density environments. In addition, the nonequilibrium state of starburst galaxies may result in short $\tau_{\rm dep}$ due to the contribution of the $\tau_{\star}/\tau_{\rm e,sf}$ term in Equation (9).

Our framework also predicts the dependence of depletion times and KSR shape on metallicity. Gas must be shielded by a certain column density in order to become cold and molecular. This column density has a corresponding number density at which such a transition occurs, as can be seen in the phase diagrams in Figure 6, that show the sharp change from the warm, transonic phase (yellow) to the cold, supersonic phase (blue) at $n\sim 10{\rm~{}cm^{-3}}$. At lower metallicities, both the characteristic number density and column density of the transition increase, leading to the decrease of $f_{\rm H_{2}}$, $f_{\rm sf}$ and $\tau_{\star}$. Thus, the overall gas depletion time, $\tau_{\rm dep}$, increases. The mechanism of $f_{\rm sf,H_{2}}$ regulation by feedback-driven turbulence discussed in Section 4.4, however, should operate regardless of the gas metallicity. Thus, we expect the KSR for molecular gas to remain linear at low metallicity. Although we also expect its normalization to be approximately independent of gas metallicity, a weak dependence is possible due to the changing properties of star-forming regions with metallicity.

The higher characteristic density of the transition at lower metallicity also results in the shift of the turnover in the KSR for total gas to higher surface densities. This shift is qualitatively similar to that predicted by the Gnedin.Kravtsov.2011 models, where star formation is tied to molecular gas. We have confirmed this explicitly by resimulating our model galaxy at a lower metallicity and will present these results in a forthcoming paper.

5.2. Implications for galaxy formation simulations

Recent studies show that in simulations with strong feedback, the galaxy-scale star formation rate — and hence the global depletion time — is insensitive to the local star formation efficiency (e.g., Dobbs.etal.2011; Agertz.etal.2013; Hopkins.etal.2013; Hopkins.etal.2017; Agertz.Kravtsov.2015; Benincasa.etal.2016; Orr.etal.2017). This behavior is thought to be due to “self-regulation” of star formation by feedback (e.g., Dobbs.etal.2011).

Our framework naturally explains the physical mechanism of such self-regulation. Equation (9) accounts for feedback via the term proportional to the mass-loading factor, $\xi$. When feedback is sufficiently strong so that this term dominates, the depletion time is given by $\tau_{\rm dep}\approx\xi\tau_{+}$, where $\tau_{+}$ is the timescale of the star-forming gas supply from the ISM. Physically, $\tau_{+}$ is not related to the local star formation efficiency, while $\xi$ is the ratio of the star-forming gas depletion time, $\tau_{\star}=\langle\epsilon_{\rm ff}/t_{\rm ff}\rangle^{-1}_{\rm sf}$, and the timescale with which feedback disperses star-forming gas, $\tau_{\rm-,fb}$ (see Section 2). Timescales $\tau_{\star}$ and $\tau_{\rm-,fb}$ are proportional to each other, because the rate of gas removal from star-forming regions is proportional to the rate of energy and momentum injection by feedback, which, in turn, is proportional to the star formation rate: ${M}_{\rm sf}/\tau_{\rm-,fb}=\xi\dot{M}_{\star}=\xi{M}_{\rm sf}/\tau_{\star}$. Thus, in this regime, the global star formation rate, $\dot{M}_{\star}={M}_{\rm g}/\tau_{\rm dep}\approx{M}_{\rm g}/(\xi\tau_{+})$, does not depend on the local depletion time and thus on local star formation efficiency. At the same time, $\dot{M}_{\star}$ is inversely proportional to the strength of feedback, $\xi$, as observed in simulations with self-regulation (e.g., Hopkins.etal.2017; Orr.etal.2017).

Note that for a given strength of feedback, the regime of self-regulation exists only for a certain range of $\epsilon_{\rm ff}$ and $\tau_{\star}$. For sufficiently low $\epsilon_{\rm ff}$ (long $\tau_{\star})$, the second term in the sum in Equation (9), i.e. $\tau_{\star}$, will dominate, and the global depletion time will scale with the local star formation efficiency as $\tau_{\rm dep}\sim\tau_{\star}\propto\epsilon_{\rm ff}^{-1}$. Indeed, we have checked that when the feedback parameters and star formation threshold in our simulations are fixed, the global depletion time is insensitive to the variation of local efficiency for $\epsilon_{\rm ff}>1\%$, but scales as $\tau_{\rm dep}\propto\epsilon_{\rm ff}^{-1}$ when $\epsilon_{\rm ff}<1\%$.

The critical value of $\epsilon_{\rm ff}$, above which $\tau_{\rm dep}$ becomes insensitive to $\epsilon_{\rm ff}$, depends on the relative contribution of the $\tau_{\star}$ term in the sum of Equation (9), and thus depends on the factors that control the first term in the sum, in particular the feedback strength, $\xi$. Dependence of the critical $\epsilon_{\rm ff}$ on feedback strength explains the results of Agertz.etal.2013 and Agertz.Kravtsov.2015, who found that the normalization of the KSR, i.e. $\tau_{\rm dep}^{-1}$, scales with $\epsilon_{\rm ff}$ in simulations without feedback but is similar for $\epsilon_{\rm ff}=1\%$ and $10\%$ when feedback is strong.

The weak sensitivity of the SFR to the value of $\epsilon_{\rm ff}$ in the regime of strong feedback does not mean that this value is not important. For example, the mass fraction of star-forming gas in this regime scales inversely with $\epsilon_{\rm ff}$: $f_{\rm sf}\equiv{M}_{\rm sf}/{M}_{\rm g}=\tau_{\star}/\tau_{\rm dep}\propto\epsilon_{\rm ff}^{-1}$. Thus, an incorrect $\epsilon_{\rm ff}$ will result in an incorrect mass fraction and an incorrect spatial distribution of star-forming gas. A similar conclusion was drawn by Hopkins.etal.2013b, who suggested that the amount of dense gas probed by HCN decreases with increasing star formation efficiency and feedback strength.

The mass fraction of star-forming gas and its distribution is also a nontrivial function of the feedback strength and the specific mix of processes that define it (e.g., Hopkins.etal.2013b; Butler.etal.2017). Stronger feedback shortens the time that gas spends in the star-forming state, $t_{\rm sf}$, and, therefore, decreases the mass fraction of star-forming gas: $f_{\rm sf}\equiv{M}_{\rm sf}/{M}_{\rm g}\sim t_{\rm sf}/(t_{\rm nsf}+t_{\rm sf})$. Besides, strong feedback shapes the overall gas PDF and may increase timescales $t_{\rm nsf}$, which also results in a decrease of $f_{\rm sf}$. This happens, for example, when a gaseous disk is stabilized by a pressure corresponding to a slowly dissipating energy directly injected by feedback (e.g., Agertz.etal.2013; Benincasa.etal.2016).

Finally, $t_{\rm nsf}$, $t_{\rm sf}$ and $f_{\rm sf}$ depend on the overall adopted definition of star-forming gas. Usually, the boundary between star-forming and non-star-forming gas is established by a set of thresholds in various gas properties, such as gas density, molecular fraction, virial parameter $\alpha_{\rm vir}$, etc. Our results indicate that these thresholds should be carefully chosen, because if thresholds grossly misidentify the star-forming gas, the values of $\epsilon_{\rm ff}$ and feedback strength may need to compensate for the wrong choice to get reasonable $t_{\rm sf}$ and $f_{\rm sf}$ values. It is not yet clear whether such compensation is possible in general, but, in any case, wrong thresholds may drive the $\epsilon_{\rm ff}$ and feedback parameters to incorrect values.

In addition, the choice of the thresholds defining star-forming gas indirectly affects the efficiency of feedback. When more gas is designated as star-forming for a given fixed star formation rate, the feedback mass-loading factor, $\xi$, will have to be larger, i.e., feedback has to be stronger. The specific choice of the thresholds matters as well. For example, when star-forming gas is defined using a density threshold, $n_{\rm sf}$, feedback has to drive star-forming gas to lower densities, $n<n_{\rm sf}$, in order to shut down star formation. If star-forming gas is defined using a threshold in the virial parameter instead, $\alpha_{\rm vir}$ can be increased to values larger than the threshold very quickly by injecting thermal or turbulent energy without changing the gas density significantly.

Thus, the framework presented in this paper implies that to get the value of $f_{\rm sf}$ and depletion time, $\tau_{\rm dep}=\tau_{\star}/f_{\rm sf}$, and the distribution of star-forming gas correctly, the overall definition of the star-forming gas, its star formation efficiency, and the strength of stellar feedback are all important and should all be modeled carefully. Ideally, the modeling choices should be based on solid physical ground or subgrid models calibrated on higher-resolution simulations. For example, recent studies show that local $\epsilon_{\rm ff}$ can be modeled using the results of high-resolution simulations of turbulent star-forming gas (Padoan.etal.2012; Padoan.etal.2017; Federrath.2015) and a subgrid model for turbulence calibrated on turbulence simulations (e.g., Braun.etal.2014; Braun.Schmidt.2015; Semenov.etal.2016; Li.etal.2017).

Interestingly, the results presented in Section 4.4 indicate that feedback can break the direct connection between the small- and large-scale relations between gas density and SFR. Specifically, our simulations assume ${\dot{\rho}}_{\star}\propto\rho^{1.5}$ on the scales of individual cells. However, as can be seen in Figures 3 and 9, on kiloparsec scales, the relation between the surface density of molecular gas and star formation rate is almost linear. This result provides a counterexample to the arguments that in simulations, the slope of the Kennicutt-Schmidt relation on kiloparsec scales simply reflects the assumed density dependence of the local SFR on small scales (Schaye.DallaVecchia.2008; Gnedin.etal.2014). Our analysis shows that the linear relation arises because only a fraction of cold, dense gas is forming stars and its density PDF does not scale self-similarly with the large-scale surface density of molecular gas, $\Sigma_{\rm H_{2}}$. Variation in covering fraction allows $\Sigma_{\rm H_{2}}$ to vary significantly, while the average density of star-forming gas and the corresponding depletion time vary only weakly. In addition, self-regulation by feedback ties the fraction of molecular gas that is forming stars and the depletion time of star-forming gas, so that their trends with $\Sigma_{\rm H_{2}}$ are similar. These trends cancel out resulting in an almost constant $\tau_{\rm dep,H_{2}}$ and linear $\dot{\Sigma}_{\star}-\Sigma_{\rm H_{2}}$ relation.

The emergence of the large-scale KSR via nontrivial, nonlinear effects of the star formation-feedback loop motivates efforts to model such processes in high-resolution simulations rather than tuning the star formation prescription to produce a particular large-scale relation.

Our model also sheds some light on the importance of the so-called “early feedback” — a collective name for the energy and momentum injection by young stars during the first 3 Myr of the life of a stellar population, before first supernovae explode. The time lag between the formation of stellar particles and the onset of feedback is important when $\tau_{\star}$ is comparable to the lag. In this case, a substantial fraction of the gas mass can be converted into stars in the first 3 Myr before supernova feedback can limit star formation. When $\tau_{\star}$ is long, the global depletion time becomes less sensitive to the presence or absence of stellar feedback during the first 3 Myr of the star-forming stage.

6. Summary and conclusions

We present a simple and intuitive physical model that elucidates why gas depletion times in galaxies are long compared to the timescales of the processes driving the evolution of the interstellar medium. We show that the depletion time is long not because some bottleneck in the formation of star-forming regions imposes a long evolutionary timescale, but because only a small fraction of the gas mass is converted into stars during a single star-forming stage in the evolution of a gas parcel. This fraction is small due to both the short duration of the star-forming stage, as dynamical processes and stellar feedback efficiently disperse star-forming regions, and the low intrinsic star formation efficiency of dense molecular gas. A gas parcel thus must go through many cycles transitioning between non-star-forming and star-forming states before it becomes converted into stars. Hence, even though the duration of each cycle can be short, the global depletion time is long because the number of cycles is large.

Furthermore, the difference between the global and local depletion times of molecular gas in our model arises because not all of the molecular gas is actively forming stars. Non-star-forming molecular gas appears naturally if local star formation efficiency is a strong function of the virial parameter of a region, while the molecular fraction of gas is set by its ability to shield against FUV radiation and is a function of mainly gas number density and metallicity.

We illustrate our model using the results of an isolated $L_{*}$-sized disk galaxy simulation that reproduces the observed Kennicutt-Schmidt relation for both molecular and atomic gas. We discuss the predictions of our model for the dependence of the global depletion time on properties of observed galaxies and on the parameters of star formation and feedback recipes in galaxy simulations. In particular, our model explains the weak sensitivity of the global star formation rate to the assumed local star formation efficiency reported in several recent numerical studies and the physics of this “self-regulation” (see Sections 2 and 5.2).

Additional results and conclusions can be summarized as follows.

Analysis of our simulation shows that the properties of gas parcels in the ISM evolve on timescales of $\sim 10-100{\rm~{}Myr}$ under the influence of compression by the spiral arms, ISM turbulence, and SNe-driven shocks. The relative importance of these processes varies with galactocentric radius, $R$, and average surface density. At $R\gtrsim 5$ kpc, the evolution from a warm, diffuse state to a dense, cold phase is driven mainly by compression in the spiral arms, while SNe-driven shocks and large-scale ISM turbulence dominate at smaller radii.

• 2.

During an evolutionary cycle, gas spends most of the time in the non-star-forming state, $t_{\rm nsf}>t_{\rm sf}$, whereas the time spent in the star-forming state, $t_{\rm sf}$, is limited by stellar feedback and dynamical processes to $t_{\rm sf}\sim 5-15{\rm~{}Myr}$. We find that the median $t_{\rm nsf}$ varies with gas surface density as $t_{\rm nsf}\propto\Sigma_{\rm g}^{-0.5}$.

• 3.

On the resolution scale of our simulation, $40$ pc, the typical range of densities in star-forming regions is limited to $n\sim 10-10^{3}{\rm~{}cm^{-3}}$ by the interstellar FUV background and stellar feedback. We find that the resulting depletion times of star-forming gas $\tau_{\star}\sim 300-500{\rm~{}Myr}$ are consistent with the depletion times estimated for observed GMCs on these scales and exhibit only weak trends with the surface densities of total and molecular gas.

• 4.

The distributions of depletion times and lifetimes of star-forming regions in our simulations imply that a typical gas parcel has to undergo $5-100$ cycles transitioning between non-star-forming and star-forming states before converting its mass into stars.

• 5.

On kiloparsec scales, our simulation produces a nearly linear relation between the surface density of H2 and surface density of star formation rate, i.e., the molecular Kennicutt-Schmidt relation, even though a nonlinear local relation, ${\dot{\rho}}_{\star}\propto\rho^{1.5}$, is adopted for star-forming gas in simulation cells. We show that the linear relation emerges due to stellar feedback, which shapes the gas density PDF making it non-self-similar and establishes a correlation between the local depletion time of the star-forming gas and the fraction of molecular gas that is in the star-forming state.

The model for the global depletion timescale presented in this paper is a generic framework that can be applied not only to galaxies as a whole but also to individual ISM patches with sizes ranging from $\sim$ kiloparsec to a typical size of star-forming regions, $\sim 10$ pc. It can also be used to predict and interpret trends of gas depletion time with the ISM properties and redshift.

As a final comment, we note that in the context of galaxy evolution over cosmological timescales, the actual gas depletion time is often considered to be unimportant. For example, galaxies from the star-forming sequence are predicted to form stars at the rates regulated by gas accretion and gas loss in winds, because gas depletion times in such galaxies are short compared to other relevant timescales (e.g., Bouche.etal.2010; Dave.etal.2012; Lilly.etal.2013). Note, however, that dwarf galaxies and galaxies from the green valley consume gas on extremely long timescales of $\gtrsim 5-10{\rm~{}Gyr}$, and, therefore, their depletion times do affect their evolution. Moreover, at $z\gtrsim 5-6$, when the age of the universe is $\lesssim 1$ Gyr, the Gyr-long gas depletion times become comparable to the cosmological evolution timescale and will therefore play an important role in controlling the SFR during the early stages of galaxies evolution (e.g., Dekel.Mandelker.2014; Peng.Maiolino.2014). The framework for modeling gas depletion time presented in this paper thus opens a way to refine theoretical models of galaxy formation in this regime, which is particularly important in the upcoming era of the JWST.

We would like to thank the referee, Mark Krumholz, for a constructive report and many detailed comments that have improved the paper. We thank Cameron Liang, Claude-André Faucher-Giguère, Robert Feldmann, Romain Teyssier, Simon Lilly, and Marcella Carollo for useful discussions. We are also grateful to Oscar Agertz, Camille Avestruz, Benedikt Diemer, Oleg Gnedin, Hui Li, Philip Mansfield, Enrique Vázquez-Semadeni, and Tony Wong, whose valuable comments and suggestions helped to improve our paper. We also thank Adam Leroy for sharing the data from Leroy.etal.2013 that we used in Figure 3. This work was supported by a NASA ATP grant NNH12ZDA001N, NSF grant AST-1412107, and by the Kavli Institute for Cosmological Physics at the University of Chicago through grant PHY-1125897 and an endowment from the Kavli Foundation and its founder, Fred Kavli. The simulation and analyses presented in this paper have been carried out using the Midway cluster at the University of Chicago Research Computing Center, which we acknowledge for support.