Agertz \& Kravtsov 2015) have demonstrated the importance of stellar feedback in shaping the key galaxy properties, such as stellar mass, size, morphology, and star formation history. The challenge of modeling feedback in galaxy formation simulations is due to the fact that the relevant processes are not fully understood and that any implementation needs to work at the resolution limit of simulations.

In our simulation, we calculate a total number of SNe II associated with a given young stellar particle assuming the Miller-Scalo IMF (Miller \& Scalo 1979). We assume that SNe explode at a uniform rate between 3 and $40 \mathrm{Myr}$ after the stellar particle is created, with each SN releasing $10^{51} \mathrm{erg}$. We use the results of $\mathrm{SN}$ remnant simulations by Martizzi et al. (2015) to estimate the fraction of this energy that gets converted into the momentum of the expanding shell (their Equation (10)) and the fraction that is left in the form of thermal energy at the snowplow stage (their Equation (9)), when the size of the $\mathrm{SN}$ remnant reaches the resolution scale $\Delta$.

Furthermore, given that in a star-forming region the remnants of multiple SNe can collide, their coherent expansion may transform into turbulent gas motions with partial momentum cancellation. As GMCs typically have distinct filamentary structure, we expect that $\sim 1 / 3$ of the total momentum is transformed into turbulent motions. Therefore, $65 \%$ of the total radial momentum given by Equation (10) of Martizzi et al. (2015) we assign as a radial momentum to 26 neighbors of the stellar particle-host cell, while $30 \%$ is converted into the SGS turbulent energy via the term $S_{\mathrm{SN}}$ in our Equation (2). Finally, we add $5 \%$ of the initial $\mathrm{SN}$ energy to the thermal energy of the host cell.

\subsection{Galaxy Model and Numerical Resolution}

In this work we follow the evolution of a Milky-Way-like disk galaxy consisting of gaseous and stellar disks with a bulge embedded in a spherical dark matter halo. We use initial conditions for the stellar component and halo similar to those used by Agertz et al. (2013), which are also similar to the initial conditions for the disk simulations in the AGORA comparison project (Kim et al. 2014).

Specifically, the dark matter halo profile has the NFW form with a concentration of $c=10$ and circular velocity at the radius enclosing the density contrast of 200 relative to the critical density of the universe of $v_{200}=150 \mathrm{~km} \mathrm{~s}^{-1}$ (Navarro et al. 1997). For the baryonic disk we use an analytical exponential density distribution with a scale radius $r_{\mathrm{d}} \approx 3.4 \mathrm{kpc}$ and scale height $h_{\mathrm{d}}=0.1 r_{\mathrm{d}}$. The spherical stellar bulge is described by the Hernquist profile (Hernquist 1990) with a scale radius $r_{\mathrm{b}}=0.1 r_{\mathrm{d}}$ and bulge-to-total disk mass ratio $M_{\mathrm{b}} / M_{\mathrm{d}}=0.1$.

The total mass of the baryonic disk is similar to that of the Milky Way: $M_{\mathrm{d}} \approx 4.3 \times 10^{10} M_{\odot}$ with $20 \%$ of mass in the gaseous disk and the rest in the stellar component. This gas mass fraction is slightly higher than the observed current $\sim 13 \%$ in the Milky Way (Nakanishi \& Sofue 2015), but is perhaps appropriate at a somewhat earlier stage of its evolution, about 1 Gyr ago.

During the disk evolution, we adaptively refine cells according to a gas mass criterion. Specifically, we split a cell when its total gas mass exceeds $8.3 \times 10^{3} M_{\odot}$. By doing so, we reach a maximal spatial resolution of $\Delta=40 \mathrm{pc}$. This scale is similar to the typical resolution of modern cosmological zoomin simulations. In addition, this scale corresponds to the size of the largest observed GMCs. Given that the star formation prescription was calibrated on the GMC scales, we decide to stop at this relatively coarse resolution and study the properties of this recipe in the regime applicable to cosmological runs. In addition to this fiducial run, we also rerun our simulation at higher resolutions, $\Delta=20$ and $10 \mathrm{pc}$, and briefly describe these results in Sections 3.3 and 3.4. We will present more detailed comparison in a follow-up paper.

We set the initial ambient halo gas density and temperature to $10^{-6} \mathrm{~cm}^{-3}$ and $10^{4} \mathrm{~K}$, respectively. As a result, the exponential gaseous disk with a constant scale height initially is slightly off thermal equilibrium. To make the transition to the equilibrium state more gradual, we initialize the subgrid turbulent velocities in the gas uniformly at the value of $\sigma=$ $1 \mathrm{~km} \mathrm{~s}^{-1}$ to provide an additional pressure support. We checked that the resulting distribution of turbulent velocities in the evolved disk does not depend on the initial value of $\sigma$.

To further minimize transient effects related to initial disk self-adjustment, we begin the evolution of the disk with cooling, star formation, and the SGS turbulence model off. We allow the disk to settle into an equilibrium and turn on dynamical production and dissipation of the SGS turbulence after 50 Myr. Cooling is then turned on gradually between 100 and 200 Myr. After that, we evolve the disk for another $100 \mathrm{Myr}$, allowing all remaining transient effects to dissipate. At this point, at $t=300$ Myr from the beginning of the simulation, the disk is in a steady thermal and dynamical equilibrium and maintains its structure, velocity and density profiles over further evolution if the star formation remains switched off. Thus, we turn on the star formation at $t=300 \mathrm{Myr}$.

For our analysis we choose a single snapshot of the disk evolution at $t \approx 600 \mathrm{Myr}$, i.e., $300 \mathrm{Myr}$ after the star formation is turned on. This choice is a compromise between allowing the disk to evolve for an appreciable amount of time in order to settle into a steady state, on the one hand, and choosing a time significantly shorter than the global gas depletion time $\left(t_{\mathrm{dep}} \sim\right.$ $1 \mathrm{Gyr})$, on the other hand. Given that our simulation does not include explicit modeling of the cosmological mass accretion, the latter condition ensures that we examine the disk before it significantly exhausts its gas supply, which could bias our comparisons with observations.

\section{RESULTS}

The SFE in our model is predicted based on the selfconsistently evolved local density, temperature, and subgrid turbulent velocity. Thus, in this section we discuss the distribution of these quantities in our simulated disk with a special emphasis on the turbulent velocities. We also compare the star formation efficiencies and rates predicted by our turbulence-based star formation prescription to observations.

As can be seen in panels (a) and (b) of Figure 1 , by $t \approx$ $600 \mathrm{Myr}$ most of the disk volume is filled by the diffuse $(n \sim$ $\left.0.1-2 \mathrm{~cm}^{-3}\right)$ warm $\left(T \sim 10^{4} \mathrm{~K}\right)$ gas, while denser gas resides in spiral structures. The spiral arms travel around the disk, compressing the diffuse gas for certain periods of time. Their motion relative to the diffuse gas is generally supersonic, and these spiral waves are thus accompanied by shocks. As gas is being compressed, at $n \sim 5 \mathrm{~cm}^{-3}$ cooling by thermal excitation of $\mathrm{C}$ II and $\mathrm{O}$ I fine-structure lines becomes efficient and the gas rapidly cools down. This substantial cooling in spiral arms means that the pre-arm shocks are radiative and, therefore, they 