© ESO, 2015

1. Introduction

Cosmological simulations are currently reaching extraordinary levels of complexity and sophistication. Based on the so-called standard model of cosmology (SMoC), these state-of-the-art simulations are able to reproduce the main features of observed large galaxies (see, e.g., Brook et al. 2012; Zemp et al. 2012; Hopkins et al. 2013; Vogelsberger et al. 2014; Perret et al. 2014; Schaye et al. 2015), although severe disagreements with observations still persist (see, e.g., Genel et al. 2014). However, the resolution of these simulations still does not allow us to treat small-scale physics in dwarf galaxies (DGs).

Unfortunately, most of the weaknesses and problems shown by the SMoC manifest themselves in the range of masses typical for DGs (see, e.g., Kroupa 2012; Famaey & McGaugh 2012, 2013, for recent reviews). Besides full-blown cosmological simulations, it is thus necessary to run simulations of individual galaxies or small groups of galaxies. This enables us to study physical processes in greater numerical resolution and address the problems of the SMoC. The simulation of individual galaxies is therefore still extremely relevant. Also, very remarkable levels of accuracy have been reached in this field (see, e.g., Hopkins et al. 2013; Renaud et al. 2013; Roškar et al. 2013; Minchev et al. 2013; Teyssier et al. 2013, among others). The simulation of individual galaxies can also be used to parameterize feedback effects that can not be treated in cosmological simulations in detail (Creasey et al. 2013; Recchi 2014).

An obvious drawback of simulating individual galaxies is that it is not clear what initial conditions one should adopt. A common strategy is to consider a rotating, isothermal gas in equilibrium with the potential generated by a fixed distribution of stars and/or of dark matter (DM), but not with the potential generated by the gas itself (Tomisaka & Ikeuchi 1988; Silich & Tenorio-Tagle 1998; Mac-Low & Ferrara 1999; Recchi & Hensler 2013). A typical justification for neglecting gas self-gravity in DG simulations is that the mass budget of these objects should be dominated by DM. In our previous paper (Vorobyov et al. 2012), we explained in detail why this approach may produce unrealistic initial gas distributions and also considered the gas self-gravity when building equilibrium initial conditions.

As known from stellar structure studies (see, e.g., Ostriker & Mark 1968), the inclusion of gas self-gravity leads to implicit equations describing the initial equilibrium configuration, which must be solved iteratively. We did that for different halo masses, degrees of rotational support, and initial temperatures and we calculated the amount of gas prone to star formation for each of these models. Our star formation criteria were based on the Toomre instability criterium and on the Kennicutt-Schmidt empirical correlation between gas surface density and star formation rate (SFR) per unit area (see Vorobyov et al. 2012, for more details). For models that satisfy both star formation criteria, the underlying assumption was that the initial configuration is marginally stable to star formation and an external perturbing agent can trigger star formation in the corresponding parts of the galaxy.

The subsequent logical step is to use these marginally stable configurations as initial conditions for full-blown hydrodynamical simulations. For this purpose, we have developed a detailed 2+1 dimensional hydrodynamical code in cylindrical coordinates with assumed axial symmetry, which treats both the stellar and gaseous components as well as the phase transitions between them. Of course, a real galaxy may depart from such axial symmetry. However, since our initial conditions are axially symmetric and we do not explicitly consider environmental effects, the departures from axial symmetry are not expected to be significant.

This is the first of a series of papers devoted to numerically studying DGs with self-consistent initial conditions. Here, we explain in detail the basic strategy and the employed numerical hydrodynamical methods, which we extensively benchmark. Furthermore, we run a few representative models to show the overall early evolution of a typical DG. We pay particular attention to the co-evolution of stars and gas in our simulations. We follow the motion of stars by means of the so-called stellar hydrodynamical approach (Theis et al. 1992; Samland et al. 1997; Vorobyov & Theis 2006; Mitchell et al. 2013; Kulikov 2014), according to which the stellar component of a galaxy can be treated as a fluid. We develop this approach further and describe how stars and gas can be coupled by means of star formation and stellar feedback processes.

The paper is organized as follows. In Sect. 2 we describe the simulation methodology, paying specific attention to coupling the stars and gas via star formation and stellar feedback. Section 3 contains the gas cooling and heating rates employed in our modeling. The solution method of our stellar hydrodynamics equations is described in Sect. 4. The evolution of several representative models of DGs is studied in Sect. 5. Model caveats are discussed in Sect. 6. The main conclusions are summarized in Sect. 7, and we benchmark our code against a suite of test problems in the Appendix.

2. Simulation methodology

Our model DGs consist of a gas disk, a stellar component, and a fixed DM halo. The gas disk is a mixture of nine chemical elements (H, He, C, O, N, Ne, Mg, Si, and Fe), with the initial abundance of heavy elements set by the adopted initial metallicity. The stellar component consists of new stars born in the course of numerical simulations and has no pre-existing component. The DM halo has a spherically symmetric form and only contributes to the total gravitational potential of the system.

We describe the time evolution of our model galaxies using the coupled system of gas and stellar hydrodynamics equations complemented with phase transformations between stars and gas, chemical enrichment by supernova explosions, low- and intermediate-mass stars, as well as star formation feedback. We assume that individual chemical elements are dynamically well coupled with the bulk motion of the gas, which means that we only need to solve for the continuity equation for every chemical element with mass density ρ_i and do not need to solve for the dynamics of every element separately. Below, we provide a brief explanation for the key equations used to evolve our system in time.

2.1. Gas hydrodynamics

The dynamics of gas is modeled using the usual hydrodynamics equations for the gas density of ρ_g, momentum ρ_gv_g, and internal energy density ϵ complemented with the rates of mass, momentum, and energy transfer between the gas and stellar components. The model also includes the continuity equations for the mass density of every chemical element ρ_i. The corresponding equations have the following form $\begin{array}{ccc}& & \frac{\mathit{\partial }{\mathit{\rho }}_{\mathrm{g}}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{(}{\mathit{\rho }}_{\mathrm{g}}{\mathit{v}}_{\mathrm{g}}\mathrm{)}\mathrm{=}\mathrm{𝒮}\mathrm{-}\mathrm{𝒟}\mathit{,}\\ & & \frac{\mathit{\partial }{\mathit{\rho }}_{\mathrm{i}}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{(}{\mathit{\rho }}_{\mathrm{i}}{\mathit{v}}_{\mathrm{g}}\mathrm{)}\mathrm{=}{\mathrm{𝒮}}^{\mathit{i}}\mathrm{-}{\mathrm{𝒟}}^{\mathit{i}}\mathit{,}\\ & & \frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{(}{\mathit{\rho }}_{\mathrm{g}}{\mathit{v}}_{\mathrm{g}}\mathrm{)}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{(}{\mathit{\rho }}_{\mathrm{g}}{\mathit{v}}_{\mathrm{g}}\mathrm{\otimes }{\mathit{v}}_{\mathrm{g}}\mathrm{)}\mathrm{=}\mathrm{-}\mathrm{\nabla }\mathit{P}\mathrm{+}{\mathit{\rho }}_{\mathrm{g}}{\mathit{g}}_{\mathrm{gr}}\\ & & \mathrm{+}{\mathit{v}}_{\mathrm{s}}\mathrm{𝒮}\mathrm{-}{\mathit{v}}_{\mathrm{g}}\mathrm{𝒟}\mathit{,}\\ & & \frac{\mathit{\partial ϵ}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{(}\mathit{ϵ}{\mathit{v}}_{\mathrm{g}}\mathrm{)}\mathrm{=}\mathrm{-}\mathit{P}\mathrm{\nabla }\mathrm{·}{\mathit{v}}_{\mathrm{g}}\mathrm{-}\mathrm{\Lambda }\mathrm{+}\mathrm{\Gamma }\mathrm{+}{\mathrm{\Gamma }}_{\mathrm{\ast }}\mathrm{-}\mathit{ϵ}\frac{\mathrm{𝒟}}{{\mathit{\rho }}_{\mathrm{g}}}\mathrm{·}\end{array}$Here, v_g and v_s are the gas and stellar velocities, P = (γ − 1)ϵ is the gas pressure linked to the internal energy density via the ideal equation of state with the ratio of specific heats γ = 5/3, and g_gr is the gravitational acceleration due to gas, stars, and DM halo. The gas cooling Λ due to line and continuum emission and gas heating Γ due to cosmic rays and small dust grains are explained in detail in Sect. 3. The gas heating due to stellar feedback Γ_∗ is calculated as the sum of the contributions from SNII and SNIa explosions (Γ_SNII and Γ_SNIa, respectively), and also from the stellar winds Γ_sw. The last term on the r.h.s. of Eq. (4) accounts for the decrease in the internal energy due to formation of stars from the gas phase. We note that the quantity ρ_gv_g ⊗ v_g is the symmetric dyadic (a rank-two tensor) calculated according to the usual rules and ∇·(ρ_gv_g⊗v_g) is the divergence of a rank-two tensor.

The source term denotes the phase transformation of stars into the gas phase and is defined as the sum of the mass release rates (per unit volume) of all individual elements, . The mass release rate of a particular element i is defined as ${\mathrm{𝒮}}^{\mathit{i}}\mathrm{=}{\mathrm{𝒮}}_{\mathrm{SNII}}^{\mathit{i}}\mathrm{+}{\mathrm{𝒮}}_{\mathrm{SNIa}}^{\mathit{i}}\mathrm{+}{\mathrm{𝒮}}_{\mathrm{\ast }}^{\mathit{i}}\mathit{,}$(5)where ${\mathrm{𝒮}}_{\mathrm{SNII}}^{\mathit{i}}$, ${\mathrm{𝒮}}_{\mathrm{SNIa}}^{\mathit{i}}$, and ${\mathrm{𝒮}}_{\mathrm{\ast }}^{\mathit{i}}$ are the contributions due to SNII, SNIa, and also due to low- and intermediate-mass stars, respectively. The sink term describes the phase transformation of gas into stars according to the adopted star formation law. The detailed explanation of all relevant rates is provided in Sect. 2.7.

2.2. Stellar hydrodynamics

The time evolution of the stellar component is computed using the Boltzmann moment equation approach introduced by Burkert & Hensler (1988), and is further developed in application to stellar disks in Samland et al. (1997) and Vorobyov & Theis (2006). In this approach, the dynamics of stars is described by the first three moments of the collisionless Boltzmann moment equations, $\begin{array}{ccc}& & \frac{\mathit{\partial }{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\left({\mathit{\rho }}_{\mathrm{s}}{\mathit{v}}_{\mathrm{s}}\right)\mathrm{=}\mathrm{𝒟}\mathrm{-}\mathrm{𝒮}\mathit{,}\\ & & \frac{\mathit{\partial }}{\mathit{\partial t}}\left({\mathit{\rho }}_{\mathrm{s}}{\mathit{v}}_{\mathrm{s}}\right)\mathrm{+}\mathrm{\nabla }\mathrm{·}\left({\mathit{v}}_{\mathrm{s}}\mathrm{\otimes }{\mathit{\rho }}_{\mathrm{s}}{\mathit{v}}_{\mathrm{s}}\right)\mathrm{=}\mathrm{-}\mathrm{\nabla }\mathrm{·}\mathrm{\Pi }\mathrm{+}{\mathit{\rho }}_{\mathrm{s}}{\mathit{g}}_{\mathrm{gr}}\mathrm{-}{\mathit{v}}_{\mathrm{s}}\mathrm{𝒮}\mathrm{+}{\mathit{v}}_{\mathrm{g}}\mathrm{𝒟}\mathit{,} \\ & & \frac{\mathit{\partial }{\mathrm{\Pi }}_{\mathit{ij}}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{(}{\mathrm{\Pi }}_{\mathit{ij}}{\mathit{v}}_{\mathrm{s}}\mathrm{)}\mathrm{=}\mathrm{-}{\mathrm{\Pi }}_{\mathit{ik}}\mathrm{:}\mathrm{(}\mathrm{\nabla }{\mathit{v}}_{\mathrm{s}}{\mathrm{)}}_{\mathit{jk}}\mathrm{-}{\mathrm{\Pi }}_{\mathit{jk}}\mathrm{:}\mathrm{(}\mathrm{\nabla }{\mathit{v}}_{\mathrm{s}}{\mathrm{)}}_{\mathit{ik}}\\ & & \mathrm{-}{\mathit{\sigma }}_{\mathit{ij}}^{\mathrm{2}}\mathrm{𝒮}\mathrm{+}{\mathit{\sigma }}_{\mathrm{\ast }}^{\mathrm{2}}\mathrm{𝒟}\mathit{,}\end{array}$where ρ_s = m^∫fd³u is the stellar mass density, ${v}{}_{\mathrm{s}}\mathrm{=}{\mathit{\rho }}_{\mathrm{s}}^{-1}\mathit{m}{}^{\mathrm{\int }}\mathit{f}\hspace{0.17em}\mathit{u}\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{u}$ is the mean stellar velocity, and Π is the stress tensor with six components ${\mathrm{\Pi }}_{\mathit{ij}}\mathrm{\equiv }{\mathit{\rho }}_{\mathrm{s}}{\mathit{\sigma }}_{\mathit{ij}}^{\mathrm{2}}$. The stellar velocity dispersions are defined as ${\mathit{\sigma }}_{\mathit{ij}}^{\mathrm{2}}\mathrm{=}{\mathit{\rho }}^{-1}\mathit{m}{}^{\mathrm{\int }}\mathit{f}\hspace{0.17em}\mathrm{(}{\mathit{u}}_{\mathit{i}}\mathrm{-}{\mathit{v}}_{\mathit{i}}\mathrm{)}\mathrm{(}{\mathit{u}}_{\mathit{j}}\mathrm{-}{\mathit{v}}_{\mathit{j}}\mathrm{)}\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}\mathit{u}$. The function f is the distribution function of stars in the six-dimensional, position-velocity phase space f ≡ f(t,x,u) and m is the average mass of a star. The quantity ∇v_s is the gradient of a vector, the covariant expression for which is provided, e.g., in Stone & Norman (1992) and Π_ik: (∇v_s)_jk is the convolution (over index k) of two rank-two tensors. The quantity σ_∗ represents a typical stellar velocity dispersion of a newborn cluster of stars, for which we take a fiducial value of 5.0 km s^-1.

2.3. Massive stars

Massive stars that end their life with SNII explosions are main contributors to the energy budget and chemical composition of the interstellar medium. Therefore, it is important to follow their evolution as accurate as possible. We improve the description of the stellar component by solving for a separate continuity equation for the mass volume density ${\mathit{\rho }}_{\mathrm{s}}^{\mathrm{h}}$ of stars with mass > 8.0 M_⊙,$\frac{{\partial }{{\rho }}_{{s}}^{{h}}}{{\partial t}}{+}{\nabla }{·}\left({{\rho }}_{{s}}^{{h}}{{v}}_{{s}}\right){=}{{f}}_{{h}}{𝒟}{-}{{𝒮}}^{{h}}{,}$(9)where f_h is the fraction of stars with mass > 8.0 M_⊙, calculated according to the adopted IMF (see Eq. (14)), and is the rate of death of massive stars. The massive star subcomponent has the same dynamical properties as the rest of the stars, an assumption that we plan to relax in the future.

2.4. Intermediate-mass stars

Stars in the (1.0−8.0) M_⊙ mass range produce the majority of carbon and nitrogen. Moreover, binary stars in this mass range can explode as Type Ia supernovae and are thus important sources of iron. Therefore, we also solve for a separate continuity equation for the mass volume density ${\mathit{\rho }}_{\mathrm{s}}^{\mathrm{m}}$ of intermediate-mass stars, $\frac{{\partial }{{\rho }}_{{s}}^{{m}}}{{\partial t}}{+}{\nabla }{·}\left({{\rho }}_{{s}}^{{m}}{{v}}_{{s}}\right){=}{{f}}_{{m}}{𝒟}{-}{{𝒮}}^{{m}}{,}$(10)where f_m is the fraction of stars with mass 1.0 M_⊙<m< 8.0 M_⊙ calculated according to the adopted IMF and is the rate of death of intermediate-mass stars. We note that the stellar evolution of low-mass stars can be neglected during the galactic evolution. They only contribute to the gravitational potential where their density can be deduced from ρ_s (the total stellar density), ${\mathit{\rho }}_{\mathrm{s}}^{\mathrm{h}}$, and ${\mathit{\rho }}_{\mathrm{s}}^{\mathrm{m}}$.

2.5. Gravity of gas, stars, and DM halo

The gravitational potential of gas and stars is calculated self-consistently by solving for the Poisson equation, ${\mathrm{\nabla }}^{\mathrm{2}}{\mathrm{\Phi }}_{\mathrm{g}\mathrm{+}\mathrm{s}}\mathrm{=}\mathrm{4}\mathit{\pi G}\mathrm{(}{\mathit{\rho }}_{\mathrm{g}}\mathrm{+}{\mathit{\rho }}_{\mathrm{s}}\mathrm{)}\mathit{.}$(11)To accelerate the simulations, we only solve for Eq. (11) when the relative error in the currently stored gravitational potential compared to the current density distribution exceeds 10^-5 (see Stone & Norman 1992, for details).

The gravitational acceleration due to a spherically symmetric DM halo g_h is calculated as explained in Vorobyov et al. (2012). Two options are available in the code: a cored isothermal DM halo profile and a Navarro-Frenk-White (or cuspy) profile. Here we use the cored isothermal DM profile. The total gravitational acceleration is calculated as g_gr = g_h − ∇Φ_{g + s}.

2.6. Mean stellar age

To calculate the mass and energy release rates by supernovae and dying low-mass stars one needs to know the age of stars. In the stellar hydrodynamics approach, stars are a mixture of populations with different ages. Nevertheless, one may define the mean age of stars ${\overline{\mathit{\tau }}}_{\mathrm{\ast }}$, which should obey the following rules:

• 1.

  The bulk motion of stars, such as compression or rarefaction, should not affect the value of ${\overline{\mathit{\tau }}}_{\mathrm{\ast }}$. This can be achieved by solving for the Lagrangian or comoving equation for ${\overline{\mathit{\tau }}}_{\mathrm{\ast }}\mathit{,}$$\frac{\mathrm{d}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}}{\mathrm{d}\mathit{t}}\mathrm{\equiv }\frac{\mathit{\partial }{\overline{\mathit{\tau }}}_{\mathrm{\ast }}}{\mathit{\partial t}}\mathrm{+}\mathrm{(}{\mathit{v}}_{\mathrm{s}}\mathrm{·}\mathrm{\nabla }\mathrm{)}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}\mathrm{=}\mathrm{0.}$(12)

• 2.

  A newly born stellar population should rejuvenate the pre-existing one. This is achieved by updating the mean age every time step in every grid cell, using the following equation: ${\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{rjv}}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{s}}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}\mathrm{+}\mathrm{△}{\mathit{M}}_{\mathrm{s}}\hspace{0.17em}\mathrm{△}\mathit{t}}{{\mathit{M}}_{\mathrm{s}}\mathrm{+}\mathrm{△}{\mathit{M}}_{\mathrm{s}}}\mathit{,}$(13)where M_s is the stellar mass in a specific cell and △ M_s the mass of stars born in the same cell during time step △ t. Equation (13) satisfies the expected asymptotic behavior of the stellar age, which becomes small (${\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{rjv}}\mathrm{\to }\mathrm{△}\mathit{t}$) during a massive instantaneous burst (△ M_s → ∞), but remains almost unchanged (${\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{rjv}}\mathrm{\to }{\overline{\mathit{\tau }}}_{\mathrm{\ast }}$) in the quiescent phase △ M_s → 0.

• 3.

  The mean stellar age should uniformly increase with time, reflecting the overall aging of the system. This is done by adding the current time step △ t to the mean stellar age in every grid cell. Thus, in the absence of star formation and pre-existing old stellar population, the mean stellar age is equivalent to the current evolution time t.

In practice, we first solve for Eq. (12) to account for bulk motions, then we update the mean stellar age in every grid cell according to Eq. (13) to account for star formation, and finally we add the current time step △ t to the mean stellar age in every grid cell to take the aging of stellar populations into account. Since we follow the evolution of both the massive and intermediate-mass star subcomponents, their mean ages ${\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{h}}$ and ${\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}$, respectively, need to be calculated using the same procedure.

Figure 1 shows the time evolution of the stellar density ρ_s (left panel) and mean stellar age ${\overline{\mathit{\tau }}}_{\mathrm{\ast }}$ (right panel) in an idealized test problem involving the gravitational collapse of a stellar sphere with unit radius. Initially, ρ_s and ${\overline{\mathit{\tau }}}_{\mathrm{\ast }}$ are set to unity inside the sphere and to zero outside the sphere (dash-dotted lines). The stellar velocity dispersion is negligible everywhere and the SFR is set to zero. As the collapse proceeds, stellar density shows the expected behavior with a central plateau shrinking in size and growing in density. At the same time the mean stellar age increases in concordance with time, as expected in the absence of star formation. In particular, the difference between the values of ${\overline{\mathit{\tau }}}_{\mathrm{\ast }}$ inside and outside the sphere is always equal to unity, as was set initially. We emphasize here the need to use the Lagrangian equation for ${\overline{\mathit{\tau }}}_{\mathrm{\ast }}$. If the Eulerian form $\mathit{\partial }{\overline{\mathit{\tau }}}_{\mathrm{\ast }}\mathit{/}\mathit{\partial t}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}{\mathit{v}}_{\mathrm{s}}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}\mathrm{\right)}\mathrm{=}\mathrm{0}$ is used instead, the above test fails. Indeed, the dash-dot-dotted line shows the mean stellar age at t = 0.4 calculated using the Eulerian form. Evidently, the difference between the values of ${\overline{\mathit{\tau }}}_{\mathrm{\ast }}$ inside and outside the sphere is now much greater than unity, indicating a spurious aging of the stellar populations inside the sphere¹.

Finally, we note that the concept of the mean stellar age has its limitations. For instance, in the case of a constant SFR and steady galaxy configuration, ${\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{h}}$ approaches a constant value after a few tens of Myr, meaning that the supernova rates and mass return rates (see Eqs. (17) and (19)) are exclusively determined by stars whose lifetime is equal to ${\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{h}}$. This bias toward stars with single age (and mass) diminishes, when a time-varying SFR is present or stellar motions are taken into account. Nevertheless, a more sophisticated approach that takes into account a possible age spread around the mean value is desirable.

2.7. Source and sink terms

To calculate the rates at which stars return their mass (including newly synthesized elements) into the gas phase, we need to make assumptions about the initial mass function (IMF), stellar lifetimes, and nucleosynthesis rates.

We adopt the Kroupa IMF (Kroupa 2001) of the form, $\mathit{\phi }\mathrm{\left(}{\mathit{M}}_{\mathrm{\ast }}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ \mathit{A}\hspace{0.17em}{\mathit{M}}_{\mathrm{\ast }}^{-1.3}& \hspace{0.17em}\hspace{0.17em}\mathrm{if} {\mathit{M}}_{\mathrm{\ast }}\mathrm{\le }\mathrm{0.5}{\mathit{M}}_{\mathrm{\odot }}\\ \mathit{B}\hspace{0.17em}{\mathit{M}}_{\mathrm{\ast }}^{-2.3}& \hspace{0.17em}\hspace{0.17em}\mathrm{if} {\mathit{M}}_{\mathrm{\ast }}\mathit{>}\mathrm{0.5}{\mathit{M}}_{\mathrm{\odot }}\mathit{,}\end{array}$(14)where the normalization constants A and B are calculated as described in Sects. 2.7.1 and 2.7.2.

The stellar lifetimes are taken from Padovani & Matteucci (1993), $\mathit{\tau }\mathrm{\left(}{\mathit{M}}_{\mathrm{\ast }}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ \mathrm{1.2}\hspace{0.17em}{\mathit{M}}_{\mathrm{\ast }}^{-1.85}\mathrm{+}\mathrm{0.003}\mathrm{Gyr}& \hspace{0.17em}\hspace{0.17em}\mathrm{if} {\mathit{M}}_{\mathrm{\ast }}\mathrm{\ge }\mathrm{7.45}{\mathit{M}}_{\mathrm{\odot }}\\ {\mathrm{10}}^{\mathit{f}\mathrm{\left(}{\mathit{M}}_{\mathrm{\ast }}\mathrm{\right)}}\mathrm{Gyr}& \hspace{0.17em}\hspace{0.17em}\mathrm{if} {\mathit{M}}_{\mathrm{\ast }}\mathit{<}\mathrm{7.45}{\mathit{M}}_{\mathrm{\odot }}\mathit{,}\end{array}$(15)where $\mathit{f}\mathrm{\left(}{\mathit{M}}_{\mathrm{\ast }}\mathrm{\right)}\mathrm{=}\frac{\mathrm{\{}\mathrm{0.334}\mathrm{-}{\sqrt{\mathrm{1.79}\mathrm{-}\mathrm{0.2232}\mathrm{\times }\left(\mathrm{7.764}\mathrm{-}\log \mathrm{\left(}{\mathit{M}}_{\mathrm{\ast }}\mathrm{\right)}\right)}}^{\mathrm{\}}}}{\mathrm{0.1116}}\mathrm{·}$(16)In Eq. (15), we adjusted the transitional value for the stellar mass (7.45 M_⊙) to obtain a smooth time derivative dM_∗/ dτ.

2.7.1. Supernova type II rates

We assume that all stars in the $\mathrm{[}\mathrm{8.0}{\mathit{M}}_{\mathrm{\odot }}\mathrm{:}{\mathit{M}}_{\mathrm{\ast }}^{\max }\mathrm{]}$ mass range end their life cycle as type II supernovae (SNeII), where ${\mathit{M}}_{\mathrm{\ast }}^{\max }$ is the upper cut-off mass of the IMF set to 100 M_⊙ in our model. The SNII rate, i.e., the number of SNeII per unit time, is then defined as ${\mathrm{ℛ}}_{\mathrm{SNII}}\mathrm{=}\{\begin{array}{c}\\ \mathit{\phi }\mathrm{\left[}{\mathit{M}}_{\mathrm{\ast }}\mathrm{\right(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{h}}\mathrm{\left)}\mathrm{\right]}\frac{\mathrm{d}{\mathit{M}}_{\mathrm{\ast }}\mathrm{\left(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{h}}\mathrm{\right)}}{\mathrm{d}\mathit{\tau }}\\ \mathrm{if}\hspace{0.17em}\hspace{0.17em}\overline{\mathit{\tau }}\begin{array}{c}\mathit{h}\\ \mathrm{\ast }\end{array}\mathrm{\in }\mathrm{\left[}\mathit{\tau }\mathrm{\right(}\mathrm{100}{\mathit{M}}_{\mathrm{\odot }}\mathrm{\left)}\mathrm{:}\mathit{\tau }\mathrm{\right(}\mathrm{8.0}{\mathit{M}}_{\mathrm{\odot }}\mathrm{\left)}\mathrm{\right]}\mathit{,}\\ \mathrm{0}\mathrm{if}\hspace{0.17em}\hspace{0.17em}\overline{\mathit{\tau }}\begin{array}{c}\mathit{h}\\ \mathrm{\ast }\end{array}\mathrm{\ni }\mathrm{\left[}\mathit{\tau }\mathrm{\right(}\mathrm{100}{\mathit{M}}_{\mathrm{\odot }}\mathrm{\left)}\mathrm{:}\mathit{\tau }\mathrm{\right(}\mathrm{8.0}{\mathit{M}}_{\mathrm{\odot }}\mathrm{\left)}\mathrm{\right]}\mathit{,}\end{array}$(17)where ${\mathit{M}}_{\mathrm{\ast }}\mathrm{\left(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{h}}\mathrm{\right)}$ is the mass of a massive star that is about to die as SNII calculated by inverting Eq. (15). This is the SNII rate following an instantaneous burst of star formation. We therefore assume here that each star formation event generates a single stellar population (SSP). The superposition of SNII rates coming from many SSPs then approximates the SNII rate in the presence of a variable SFR. The derivative dM_∗/ dτ is calculated by inverting Eq. (15) and then differentiating it with respect to time.

The IMF in Eq. (17) has to be normalized according to the total mass of massive stars ${\mathit{M}}_{\mathrm{s}\mathit{,}\mathrm{tot}}^{\mathrm{h}}$ in each computational cell, $\underset{\mathrm{8.0}{\mathit{M}}_{\mathrm{\odot }}}{\overset{{\mathit{M}}_{\mathrm{\ast }}\mathrm{\left(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{h}}\mathrm{\right)}}{\mathrm{\int }}}\mathit{\phi }\mathrm{\left(}{\mathit{M}}_{\mathrm{\ast }}\mathrm{\right)}{\mathit{M}}_{\mathrm{\ast }}\hspace{0.17em}\mathrm{d}{\mathit{M}}_{\mathrm{\ast }}\mathrm{=}{\mathit{M}}_{\mathrm{s}\mathit{,}\mathrm{tot}}^{\mathrm{h}}\mathit{.}$(18)The upper limit in the integral is not fixed, rather it depends on the mean stellar age of massive stars in a given computational cell. This is done to take into account the fact that massive stars with lifetimes smaller than ${\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{h}}$, and hence with mass greater than ${\mathit{M}}_{\mathrm{\ast }}\mathrm{\left(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{h}}\mathrm{\right)}$, must have exploded as SNeII. Excluding those stars ensures a correct calculation of the stellar mass spectrum at any time instant². Since ${\mathit{M}}_{\mathrm{s}}^{\mathrm{h}}$ and ${\mathit{M}}_{\mathrm{\ast }}\mathrm{\left(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{h}}\mathrm{\right)}$ are time-varying quantities, the normalization constants in the IMF need to be updated every time step.

The mass return rate (per unit volume) by SNII explosions can now be defined as ${\mathrm{𝒮}}_{\mathrm{SNII}}^{\mathit{i}}\mathrm{=}{\mathrm{ℛ}}_{\mathrm{SNII}}\hspace{0.17em}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{SNII}}^{\mathit{i}}\mathrm{\left[}{\mathit{M}}_{\mathrm{\ast }}\mathrm{\right(}\overline{\mathit{\tau }}\begin{array}{c}\mathrm{h}\\ \mathrm{\ast }\end{array}{\mathrm{\left)}}^{\mathrm{\right]}}\mathit{,}$(19)where ${\mathit{\rho }}_{\mathrm{SNII}}^{\mathit{i}}\mathrm{\left[}{\mathit{M}}_{\mathrm{\ast }}\mathrm{\right(}\overline{\mathit{\tau }}\begin{array}{c}\mathrm{h}\\ \mathrm{\ast }\end{array}\mathrm{\left)}\mathrm{\right]}$ is the mass (per unit volume) released by a SNII of mass M_∗ in the form of a specific element i. The values of ρⁱ for nine considered elements (H, He, C, O, N, Ne, Mg, Si, and Fe) and four stellar metallicities (Z_∗ = 10^-4,10^-2,10^-1, and 1.0 times solar) are taken from Woosley & Weaver (1995). For stellar masses greater than 40 M_⊙, the yields are set constant to those of a 40 M_⊙ star.

Finally, we assume that each supernova releases ϵ_SN × 10⁵¹ erg in the form of thermal energy and define the rate of energy density released by SNeII as ${\mathrm{\Gamma }}_{\mathrm{SNII}}\mathrm{=}{\mathrm{10}}^{\mathrm{51}}{\mathit{ϵ}}_{\mathrm{SN}}{\mathrm{ℛ}}_{\mathrm{SNII}}\hspace{0.17em}{\mathit{V}}^{-1}\mathit{,}$(20)where V is the injection volume specified by the actual size of our numerical grid and ϵ_SN is the ejection efficiency set to 1.0 in this work (see Sect. 6 for discussion).

2.7.2. Supernova type Ia rates

We derive the rates of Type Ia supernovae following the ideas and notations laid out in Matteucci & Recchi (2001). The number of binary systems (per unit mass of the binary system M_B) that are capable of producing a SNIa explosion can be calculated as ${\mathit{N}}_{\mathrm{SNIa}}\mathrm{\left[}{\mathit{M}}_{\mathrm{2}}\mathrm{\right(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}\mathrm{\left)}\mathrm{\right]}\mathrm{=}{\mathit{A}}_{\mathrm{\ast }}\underset{{\mathit{M}}_{\mathrm{B}\mathit{,}\mathrm{1}}}{\overset{{\mathit{M}}_{\mathrm{B}\mathit{,}\mathrm{2}}}{\mathrm{\int }}}\mathit{\phi }\mathrm{\left(}{\mathit{M}}_{\mathrm{B}}\mathrm{\right)}\mathit{f}\left(\frac{{\mathit{M}}_{\mathrm{2}}\mathrm{\left(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}\mathrm{\right)}}{{\mathit{M}}_{\mathrm{B}}}\right)\frac{\mathrm{d}{\mathit{M}}_{\mathrm{B}}}{{\mathit{M}}_{\mathrm{B}}}\mathit{,}$(21)where ${\mathit{M}}_{\mathrm{2}}\mathrm{\left(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}\mathrm{\right)}$ is the mass of the secondary component, A_∗ a normalization constant and the limits of the integral are defined as $\begin{array}{ccc}{\mathit{M}}_{\mathrm{B}\mathit{,}\mathrm{1}}& \mathrm{=}& \begin{array}{c}\max \end{array}\left[\mathrm{2}{\mathit{M}}_{\mathrm{2}}\mathrm{\left(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}\mathrm{\right)}\mathit{,}{\mathit{M}}_{\mathrm{B}\mathit{,}\min }\right]\\ {\mathit{M}}_{\mathrm{B}\mathit{,}\mathrm{2}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}{\mathit{M}}_{\mathrm{B}\mathit{,}\max }\mathrm{+}{\mathit{M}}_{\mathrm{2}}\mathrm{\left(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}\mathrm{\right)}\mathit{.}\end{array}$The maximum and minimum masses of a binary system, which is capable of producing a SNIa explosion, are set to M_B,min = 3.0 M_⊙ and M_B,max = 16.0 M_⊙, respectively, and the maximum mass of the secondary is M₂ = 8.0 M_⊙. The quantity $\mathit{f}\mathrm{\left(}{\mathit{M}}_{\mathrm{2}}\mathrm{\right(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}\mathrm{)}\mathit{/}{\mathit{M}}_{\mathrm{B}}\mathrm{)}$ in Eq. (21) is the distribution function of the ${\mathit{M}}_{\mathrm{2}}\mathrm{\left(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}\mathrm{\right)}\mathit{/}{\mathit{M}}_{\mathrm{B}}$ ratio and is defined as $\mathit{f}\left(\frac{{\mathit{M}}_{\mathrm{2}}\mathrm{\left(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}\mathrm{\right)}}{{\mathit{M}}_{\mathrm{B}}}\right)\mathrm{=}{\mathrm{2}}^{\mathrm{1}\mathrm{+}\mathit{\gamma }}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\gamma }\mathrm{)}{\left(\frac{{\mathit{M}}_{\mathrm{2}}\mathrm{\left(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}\mathrm{\right)}}{{\mathit{M}}_{\mathrm{B}}}\right)}^{\mathit{\gamma }}\mathit{,}$(24)where γ is set to 0.5. For this value of γ, the normalization constant that better reproduces the observed [α/Fe] ratios in dwarf galaxies is A_∗ = 0.06 (Recchi et al. 2009), therefore we adopt this value.

By analogy to the massive star subcomponent, the IMF in Eq. (21) has to be normalized according to the total mass of intermediate-mass stars ${\mathit{M}}_{\mathrm{s}\mathit{,}\mathrm{tot}}^{\mathrm{m}}$ in each computational cell, $\underset{\mathrm{1.0}{\mathit{M}}_{\mathrm{\odot }}}{\overset{{\mathit{M}}_{\mathrm{\ast }}\mathrm{\left(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}\mathrm{\right)}}{\mathrm{\int }}}\mathit{\phi }\mathrm{\left(}{\mathit{M}}_{\mathrm{\ast }}\mathrm{\right)}{\mathit{M}}_{\mathrm{\ast }}\hspace{0.17em}\mathrm{d}{\mathit{M}}_{\mathrm{\ast }}\mathrm{=}{\mathit{M}}_{\mathrm{s}\mathit{,}\mathrm{tot}}^{\mathrm{m}}\mathit{.}$(25)The SNIa rate, i.e., the number of SNIa explosions per unit time, can now be calculated by analogy to Eq. (17) as ${\mathrm{ℛ}}_{\mathrm{SNIa}}\mathrm{=}\{\begin{array}{c}\\ {\mathit{N}}_{\mathrm{SNIa}}\left[{\mathit{M}}_{\mathrm{2}}\mathrm{\left(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}\mathrm{\right)}\mathrm{)}\right]\frac{\mathrm{d}{\mathit{M}}_{\mathrm{2}}\mathrm{\left(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}\mathrm{\right)}\mathrm{)}}{\mathrm{d}\mathit{\tau }}\\ \mathrm{if}\hspace{0.17em}\hspace{0.17em}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}\mathrm{\in }\mathrm{\left[}\mathit{\tau }\mathrm{\right(}\mathrm{8.0}{\mathit{M}}_{\mathrm{\odot }}\mathrm{\left)}\mathrm{:}\mathit{\tau }\mathrm{\right(}{\mathit{M}}_{\mathrm{2}\mathit{,}\min }\mathrm{\left)}\mathrm{\right]}\mathit{,}\\ \mathrm{0}\mathrm{if}\hspace{0.17em}\hspace{0.17em}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}\mathrm{\ni }\mathrm{\left[}\mathit{\tau }\mathrm{\right(}\mathrm{8.0}{\mathit{M}}_{\mathrm{\odot }}\mathrm{\left)}\mathrm{:}\mathit{\tau }\mathrm{\right(}{\mathit{M}}_{\mathrm{2}\mathit{,}\min }\mathrm{\left)}\mathrm{\right]}\mathit{,}\end{array}$(26)where M_2,min = 1.0M_⊙ is the minimum mass of the secondary.

The mass return rate (per unit volume) by SNIa explosions can now be expressed as ${\mathrm{𝒮}}_{\mathrm{SNIa}}^{\mathit{i}}\mathrm{=}{\mathrm{ℛ}}_{\mathrm{SNIa}}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{SNIa}}^{\mathit{i}}\mathrm{\left[}{\mathit{M}}_{\mathrm{\ast }}\mathrm{\right(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}\mathrm{\left)}\mathrm{\right]}\mathit{,}$(27)where ${\mathit{\rho }}_{\mathrm{SNIa}}^{\mathit{i}}\mathrm{\left[}{\mathit{M}}_{\mathrm{\ast }}\mathrm{\right(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}\mathrm{\left)}\mathrm{\right]}$ is the mass (per unit volume) released by a SNIa of mass M_∗ in the form of a specific element i. The values of ρⁱ for six elements (C, O, N, Mg, Si, and Fe) are taken from Nomoto et al. (1984). We neglect the metallicity dependence of the SNIa yields. The energy release rate by SNIa (Γ_SNIa) is defined in analogy to Eq. (20), with ℛ_SNII being substituted by ℛ_SNIa.

2.7.3. Stellar winds and the mass return by intermediate-mass stars

The luminosity produced by stellar winds is taken from results of the Starburst99 software package (Leitherer et al. 1999; Vazquez & Leitherer 2005). We used Padova AGB stellar tracks at different metallicities (Z_∗ = 0.0004, 0.004, 0.008, 0.02, 0.05) and used this set of models as a library to obtain wind luminosities for each SSP, depending on its mass and average metallicity. The energy transfer efficiency of the wind power is set to 100%.

The mass-return rate (per unit volume) by winds of stars in the (1.0−8.0) M_⊙ mass range is calculated as ${\mathrm{𝒮}}_{\mathrm{\ast }}^{\mathit{i}}\mathrm{=}{\mathrm{ℛ}}_{\mathrm{\ast }}{\mathit{\rho }}_{\mathrm{\ast }}^{\mathit{i}}\left[{\mathit{M}}_{\mathrm{\ast }}\mathrm{\left(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}\mathrm{\right)}\right]\mathit{,}$(28)where the number of stars with mass 1.0 <M_∗< 8.0 M_⊙ dying per unit time is expressed as ${\mathrm{ℛ}}_{\mathrm{\ast }}\mathrm{=}\{\begin{array}{c}\\ \mathit{\phi }\mathrm{\left[}{\mathit{M}}_{\mathrm{\ast }}\mathrm{\right(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}\mathrm{\left)}\mathrm{\right]}\frac{\mathrm{d}{\mathit{M}}_{\mathrm{\ast }}\mathrm{\left(}{\overline{\mathit{\tau }}}_{\mathrm{\ast }}^{\mathrm{m}}\mathrm{\right)}}{\mathrm{d}\mathit{\tau }}\\ \mathrm{if}\hspace{0.17em}\hspace{0.17em}\mathit{t}\mathrm{\in }\mathrm{\left[}\mathit{\tau }\mathrm{\right(}\mathrm{8.0}{\mathit{M}}_{\mathrm{\odot }}\mathrm{\left)}\mathrm{:}\mathit{\tau }\mathrm{\right(}\mathrm{1.0}{\mathit{M}}_{\mathrm{\odot }}\mathrm{\left)}\mathrm{\right]}\mathit{,}\\ \mathrm{0}\mathrm{if}\hspace{0.17em}\hspace{0.17em}\mathit{t}\mathrm{\ni }\mathrm{\left[}\mathit{\tau }\mathrm{\right(}\mathrm{8.0}{\mathit{M}}_{\mathrm{\odot }}\mathrm{\left)}\mathrm{:}\mathit{\tau }\mathrm{\right(}\mathrm{1.0}{\mathit{M}}_{\mathrm{\odot }}\mathrm{\left)}\mathrm{\right]}\mathit{,}\end{array}$(29)and the mass (per unit volume) returned by dying stars ${\mathit{\rho }}_{\mathrm{\ast }}^{\mathit{i}}$ in the form of a specific element i includes both the true nucleosynthesis yields of C, O, and N and the pre-existed contribution (see Recchi et al. 2001). The nucleosynthesis yields are taken from van den Hoek & Groenewegen (1997).

2.7.4. Star formation rate

The phase transformation of gas into stars is controlled by the following SFR per unit volume: $\begin{array}{ccc}\mathrm{𝒟}\mathrm{=}\{\begin{array}{c}\\ \\ \end{array}& & \end{array}$(30)where n_g and T are the gas number density and temperature, respectively, ∇·v_g is the gas velocity divergence, and ξ_cpw is the cold plus warm gas fraction (see Sect. 3.4 for detail). The exponent 1.5 in Eq. (30) is expected if the SFR is proportional to the ratio of the local gas volume density to the free-fall time, and the free-fall and cooling times are calculated as ${\mathit{\tau }}_{\mathrm{ff}}\mathrm{=}\sqrt{\mathrm{3}\mathit{\pi }\mathit{/}\mathrm{\left(}\mathrm{32}\mathit{G}{\mathit{\rho }}_{\mathrm{g}}\mathrm{\right)}}$ and τ_c = ϵ/ Λ, respectively³. The term in brackets applies the adopted temperature dependence by quickly turning off star formation at temperatures greater than T₀ = 10³ K. For consistency with other studies (e.g., Springel & Hernquist 2003), the critical gas number density above which star formation is allowed is set to 1.0 cm^-3 (see discussion in Sect. 6). Finally, the normalization constant c_∗ is set to 0.05 (Stinson et al. 2006).

2.8. Stellar metallicity

The heavy element yields by SNeII and also by low- and intermediate-mass stars are metallicity dependent. We calculate the stellar metallicity Z_∗ using a procedure similar to that applied to the mean stellar age in Sect. 2.6. To account for bulk motions, we first solve for the Lagrangian equation for Z_∗ written in the following form: $\frac{\mathrm{d}{\mathit{Z}}_{\mathrm{\ast }}}{\mathrm{d}\mathit{t}}\mathrm{\equiv }\frac{\mathit{\partial }{\mathit{Z}}_{\mathrm{\ast }}}{\mathit{\partial t}}\mathrm{+}\mathrm{(}{\mathit{v}}_{\mathrm{s}}\mathrm{·}\mathrm{\nabla }\mathrm{)}{\mathit{Z}}_{\mathrm{\ast }}\mathrm{=}\mathrm{0.}$(31)Then, we update the stellar metallicity in every grid cell to account for the production of metals due to star formation using the following equation: ${\mathit{Z}}_{\mathrm{\ast }}^{\mathrm{new}}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{s}}{\mathit{Z}}_{\mathrm{\ast }}\mathrm{+}\mathrm{△}{\mathit{M}}_{\mathrm{s}}\hspace{0.17em}{\mathit{Z}}_{\mathrm{g}}}{{\mathit{M}}_{\mathrm{s}}\mathrm{+}\mathrm{△}{\mathit{M}}_{\mathrm{s}}}\mathit{,}$(32)where Z_g is the current metallicity of gas from which stars form. Equation (32) satisfies the expected asymptotic behavior of the stellar metallicity, which approaches that of the gas (${\mathit{Z}}_{\mathrm{\ast }}^{\mathrm{new}}\mathrm{\to }{\mathit{Z}}_{\mathrm{g}}$) during a massive instantaneous burst (△ M_s → ∞), but remains almost unchanged (${\mathit{Z}}_{\mathrm{\ast }}^{\mathrm{new}}\mathrm{\to }{\mathit{Z}}_{\mathrm{\ast }}$) in the quiescent phase △ M_s → 0. Step 3 from Sect. 2.6 does not need to be applied to the stellar metallicity because this step takes the overall aging of stellar populations into account.

3. Cooling and heating rates

We assume that the gas is in collisional ionization equilibrium and calculate the gas cooling rates using the cooling functions presented in Böhringer & Hensler (1989) and Dalgarno & McCray (1972).

3.1. Low-temperature cooling

For low-temperature cooling T< 10⁴ K, we adopted the cooling rates as described in Dalgarno & McCray (1972). The chemical elements that are taken into account are H, O, C, N, Si, and Fe. For high fractional ionization, cooling is mostly provided by collisions of neutral and singly-ionized ions with thermal electrons. We assumed that most of oxygen and nitrogen is in neutral form. The following equations describe the cooling efficiency L_e(X_i) (in erg cm³ s^-1) due to collisions of element X_i of number density n_i with thermal electrons of number density n_e: $\begin{array}{ccc}& & {\mathit{L}}_{\mathrm{e}}^{\mathrm{low}}\mathrm{\left(}{\mathrm{C}}^{\mathrm{+}}\mathrm{\right)}\mathrm{=}{\mathit{T}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{(}\mathrm{7.9}\mathrm{\times }{\mathrm{10}}^{-20}{\mathrm{e}}^{\mathrm{-}\mathrm{92}\mathit{/}\mathit{T}}\mathrm{+}\mathrm{3.0}\mathrm{\times }{\mathrm{10}}^{-17}{\mathrm{e}}^{\mathrm{-}\mathrm{61900}\mathit{/}\mathit{T}}\mathrm{)}\mathit{,} \\ & & {\mathit{L}}_{\mathrm{e}}^{\mathrm{low}}\mathrm{\left(}{\mathrm{Si}}^{\mathrm{+}}\mathrm{\right)}\mathrm{=}{\mathit{T}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{(}\mathrm{1.9}\mathrm{\times }{\mathrm{10}}^{-18}{\mathrm{e}}^{\mathrm{-}\mathrm{413}\mathit{/}\mathit{T}}\mathrm{+}\mathrm{3.0}\mathrm{\times }{\mathrm{10}}^{-17}{\mathrm{e}}^{\mathrm{-}\mathrm{63600}\mathit{/}\mathit{T}}\mathrm{)}\mathit{,} \\ & & {\mathit{L}}_{\mathrm{e}}^{\mathrm{low}}\mathrm{\left(}{\mathrm{Fe}}^{\mathrm{+}}\mathrm{\right)}\mathrm{=}{\mathit{T}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{(}\mathrm{4.8}\mathrm{\times }{\mathrm{10}}^{-18}{\mathrm{e}}^{\mathrm{-}\mathrm{2694}\mathit{/}\mathit{T}}\\ & & \mathrm{+}\mathrm{1.1}\mathrm{\times }{\mathrm{10}}^{-18}\mathrm{[}{\mathrm{e}}^{\mathrm{-}\mathrm{554}\mathit{/}\mathit{T}}\mathrm{+}\mathrm{1.3}{\mathrm{e}}^{\mathrm{-}\mathrm{961}\mathit{/}\mathit{T}}\mathrm{]}\\ & & \mathrm{+}\mathrm{7.8}\mathrm{\times }{\mathrm{10}}^{-18}{\mathrm{e}}^{\mathrm{-}\mathrm{3496}\mathit{/}\mathit{T}}\mathrm{)}\mathit{,}\\ & & {\mathit{L}}_{\mathrm{e}}^{\mathrm{low}}\mathrm{\left(}\mathrm{O}\mathrm{\right)}\mathrm{=}\mathrm{1.74}\mathrm{\times }{\mathrm{10}}^{-24}{\mathit{T}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{[}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{7.6}{\mathit{T}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{)}{\mathrm{e}}^{\mathrm{-}\mathrm{228}\mathit{/}\mathit{T}}\\ & & \mathrm{+}\mathrm{0.38}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{7.7}{\mathit{T}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{)}{\mathrm{e}}^{\mathrm{-}\mathrm{326}\mathit{/}\mathit{T}}\mathrm{]}\\ & & \mathrm{+}\mathrm{9.4}\mathrm{\times }{\mathrm{10}}^{-23}{\mathit{T}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathrm{e}}^{\mathrm{-}\mathrm{22700}\mathit{/}\mathit{T}}\mathit{,}\\ & & {\mathit{L}}_{\mathrm{e}}^{\mathrm{low}}\mathrm{\left(}\mathrm{N}\mathrm{\right)}\mathrm{=}\mathrm{8.2}\mathrm{\times }{\mathrm{10}}^{-22}{\mathit{T}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{2.7}\mathrm{\times }{\mathrm{10}}^{-9}{\mathit{T}}^{\mathrm{2}}\mathrm{)}{\mathrm{e}}^{\mathrm{-}\mathrm{27700}\mathit{/}\mathit{T}}\mathit{.}\end{array}$The cooling rate per unit volume Λ^low (erg cm^-3 s^-1) due to collisions of element X_i with thermal electrons is then calculated as ${\mathit{L}}_{\mathrm{e}}^{\mathrm{low}}\mathrm{\left(}{\mathit{X}}_{\mathit{i}}\mathrm{\right)}\hspace{0.17em}{\mathit{n}}_{\mathrm{i}}\hspace{0.17em}{\mathit{n}}_{\mathrm{e}}$. We also take cooling due to collisions of hydrogen atoms with thermal electrons, as tabulated in Table 2 of Dalgarno & McCray (1972) into account.

For low fractional ionization, collisions of C⁺, O, Si⁺, and Fe⁺ with neutral hydrogen may become important contributors to the total cooling budget. The following equations describe the corresponding cooling efficiency for Si⁺ and Fe⁺:$\begin{array}{ccc}{\mathit{L}}_{\mathrm{H}}^{\mathrm{low}}\mathrm{\left(}{\mathrm{Si}}^{\mathrm{+}}\mathrm{\right)}& \mathrm{=}& \mathrm{7.4}\mathrm{\times }{\mathrm{10}}^{-23}{\mathrm{e}}^{\mathrm{-}\mathrm{413}\mathit{/}\mathit{T}}\mathit{,}\\ {\mathit{L}}_{\mathrm{H}}^{\mathrm{low}}\mathrm{\left(}{\mathrm{Fe}}^{\mathrm{+}}\mathrm{\right)}& \mathrm{=}& \mathrm{1.1}\mathrm{\times }{\mathrm{10}}^{-22}\mathrm{(}{\mathrm{e}}^{\mathrm{-}\mathrm{554}\mathit{/}\mathit{T}}\mathrm{+}\mathrm{1.4}{\mathrm{e}}^{\mathrm{-}\mathrm{961}\mathit{/}\mathit{T}}\mathrm{)}\mathit{.}\end{array}$The cooling rate per unit volume due to collisions of element X_i with hydrogen atoms of number density n_H is then calculated as ${\mathit{L}}_{\mathrm{H}}^{\mathrm{low}}\mathrm{\left(}{\mathit{X}}_{\mathit{i}}\mathrm{\right)}\hspace{0.17em}{\mathit{n}}_{\mathit{i}}\hspace{0.17em}{\mathit{n}}_{\mathrm{H}}$. The cooling efficiencies due to collisions of O and C⁺ with neutral hydrogen are taken from Table 4 of Dalgarno & McCray (1972).

To summarize, the total cooling rate at the low-temperature regime is calculated as $\begin{array}{ccc}{\mathrm{\Lambda }}^{\mathrm{low}}\left[\frac{\mathrm{erg}}{\mathrm{c}{\mathrm{m}}^{\mathrm{3}}\mathrm{s}}\right]\mathrm{=}\sum _{\mathit{i}}{\mathit{n}}_{\mathit{i}}{\mathit{n}}_{\mathrm{H}}{\mathit{\xi }}_{\mathrm{cpw}}\mathrm{[}{\mathit{L}}_{\mathrm{e}}^{\mathrm{low}}\mathrm{\left(}{\mathit{X}}_{\mathit{i}}\mathrm{\right)}\hspace{0.17em}{\mathit{f}}_{\mathrm{ion}}\mathrm{+}{\mathit{L}}_{\mathrm{H}}^{\mathrm{low}}\mathrm{\left(}{\mathit{X}}_{\mathit{i}}\mathrm{\right)}{\hspace{0.17em}}^{\mathrm{]}}\mathit{,}& & \end{array}$(40)where summation is done over six elements: H, O, C, N, Si, and Fe, where applicable. Here, ξ_cpw is the cold plus warm gas fraction, and f_ion = n_e/n_H is the ionization fraction, the value of which for a given gas temperature T_g is taken from Table 2 of Schure et al. (2009) for the collisionally ionized gas in thermal equilibrium for T_g> 10^3.8 K. At lower temperatures, f_ion is set to a constant value of 0.01.

Most of the considered elements (C, Si, Fe, and O) have relatively low critical densities n_crit at which spontaneous and collisional de-excitations become comparable for some transition levels. Cooling due to these elements at n_H ≳ n_crit is proportional to n_in_crit rather than to n_in_H. We take this into account by multiplying Eq. (40) with n_crit/ (n_crit + n_H) and taking n_crit from Table 8 in Hollenbach & McKee (1989). At a low numerical resolution, the critical densities may never be reached.

3.2. High-temperature cooling

For the solar-metallicity plasma with temperature T ≥ 10⁴ K, we adopted the cooling rates of Böhringer & Hensler (1989), which include the line emission for most abundant elements and continuum emission due to free-free, two-photon, and recombination radiation. We have considered nine chemical elements (H, He, C, N, O, Ne, Mg, Si, and Fe), the total (line plus continuum) cooling efficiencies L^high(X_i) for which are plotted in Fig. 2 of Böhringer & Hensler (1989). These cooling efficiencies are normalized to n_en_H and hence need to be multiplied by the ionization fraction f_ion = n_e/n_H to retrieve the values that can later be used in our hydro code. The total cooling rates in the high-temperature regime is then calculated as ${\mathrm{\Lambda }}^{\mathrm{high}}\left[\frac{\mathrm{erg}}{\mathrm{c}{\mathrm{m}}^{\mathrm{3}}\mathrm{s}}\right]\mathrm{=}\sum _{\mathit{i}}{\mathit{L}}^{\mathrm{high}}\mathrm{\left(}{\mathit{X}}_{\mathit{i}}\mathrm{\right)}\hspace{0.17em}{\mathit{f}}_{\mathrm{ion}}\hspace{0.17em}\mathrm{(}{\mathit{n}}_{\mathrm{H}}{\mathit{\xi }}_{\mathrm{cpw}}{\mathrm{)}}^{\mathrm{2}}\frac{{\mathit{n}}_{\mathit{i}}\hspace{0.17em}}{{\mathit{n}}_{\mathit{i}}^{\mathrm{sol}}}\mathit{,}$(41)where n_i is the current number density of element i, ${\mathit{n}}_{\mathit{i}}^{\mathrm{sol}}$ is its number density at the solar abundance, and summation is performed over nine elements: H, He, C, N, O, Ne, Mg, Si, and Fe.

Figure 2 shows the cooling efficiency laid out in the following form, $\sum _{\mathit{i}}\frac{{\mathit{n}}_{\mathit{i}}}{{\mathit{n}}_{\mathit{H}}}\mathrm{\left[}{\mathit{f}}_{\mathrm{ion}}{\mathit{L}}_{\mathrm{e}}^{\mathrm{low}}\mathrm{\right(}{\mathit{X}}_{\mathit{i}}\mathrm{)}\mathrm{+}{\mathit{L}}_{\mathrm{H}}^{\mathrm{low}}\mathrm{(}{\mathit{X}}_{\mathit{i}}{\mathrm{\left)}}^{\mathrm{\right]}}\mathrm{+}{\mathit{f}}_{\mathrm{ion}}{\mathit{L}}^{\mathrm{high}}\mathrm{\left(}{\mathit{X}}_{\mathit{i}}\mathrm{\right)}\mathit{.}$(42)In particular, the left panel has f_ion = 1, but with varying metallicity, and the right panel has the solar metallicity, but with varying f_ion. A strong dependence on the metallicity and ionization fraction for f_ion ≳ 10^-3 is evident in the figure.

3.3. Gas temperature

To use the cooling rates described above, one needs to know the gas temperature T. The latter is calculated from the ideal equation of state P = ρ_gℛT_g/μ, where ℛ is the universal gas constant, using the following expression for the mean molecular weight: $\mathit{\mu }\mathrm{=}\frac{\sum _{\mathit{i}}{\mathit{\rho }}_{\mathit{i}}\mathrm{+}{\mathit{m}}_{\mathrm{e}}{\mathit{f}}_{\mathrm{ion}}{\mathit{n}}_{\mathrm{H}}}{\sum _{\mathit{i}}{\mathit{n}}_{\mathit{i}}\mathrm{+}{\mathit{f}}_{\mathrm{ion}}{\mathit{n}}_{\mathrm{H}}}\mathit{,}$(43)where m_e is the electron mass and the summation is performed over all nine chemical elements considered in this study.

3.4. Avoiding the overcooling of SN ejecta

When the supernova energy is released in the form of the thermal energy, there may exist the risk of overestimating the gas cooling, which may lead to radiating away most of the released thermal energy. This occurs because the cooling rate Λ is proportional to the square of gas number density, which, for a single-phase medium, may be significantly higher than what is expected in the case of a multiphase medium for the hot supernova ejecta.

To avoid the spurious gas overcooling, we make use of the following procedure. We keep track of the hot SN ejecta by solving for the following equation for the hot gas density ρ_g,h: $\frac{{\partial }{{\rho }}_{{g}{,}{h}}}{{\partial t}}{+}{\nabla }{·}\left({{\rho }}_{{g}{,}{h}}{{v}}_{{g}}\right){=}{𝒮}{-}\frac{{{\rho }}_{{g}{,}{h}}}{{{\tau }}_{{c}{,}{h}}}{·}$(44)The second term on the r.h.s. accounts for the phase transformation of the hot gas into the cold plus warm (CpW) component, where τ_c,h = ϵ/ Λ_h denotes the cooling time of the hot component and the cooling rate Λ_h is calculated assuming the characteristic cooling efficiency of 5 × 10^-23 ergs cm³ s^-1 of the hot component with the typical temperature of 10⁸ K (see Fig. 2 and Sect. 6)⁴.

Equation (44) is solved along with the system of Eqs. (1)−(4) describing the evolution of the gas disk, and the CpW gas fraction is calculated as ξ_cpw = (ρ_g − ρ_g,h) /ρ_g. We then modify the cooling rates and the SFR by multiplying the total gas density ρ_g and the hydrogen number density n_H in Eqs. (40) and (41) and also in Eq. (30) with ξ_cpw to retrieve the density of the CpW component. This procedure scales down the cooling rates and the SF rate according to the mass of the CpW component in each computational cell and helps to avoid overcooling and unphysically high SF rates.

3.5. Heating processes

3.5.1. Cosmic ray heating

The cosmic ray heating efficiency is the product of the cosmic ray ionization rate ξ_cr and the energy input per ionization △ Q_cr. We adopt ξ_cr = 4 × 10^-17 s^-1 based on the work of Van Dishoeck & Black (1986). The value of △ Q_cr is approximately 20 eV. The resulting heating rate due to ionization of hydrogen and helium atoms is ${\mathrm{\Gamma }}_{\mathrm{cr}}\mathrm{=}\mathrm{1.34}\mathrm{\times }{\mathrm{10}}^{-27}{\mathit{G}}_{\mathrm{0}}\left(\mathrm{0.48}{\mathit{n}}_{\mathrm{H}}\mathrm{+}\mathrm{0.5}{\mathit{n}}_{\mathrm{He}}\right)\left[\frac{\mathrm{ergs}}{{\mathrm{cm}}^{\mathrm{3}}\mathrm{s}}\right]\mathit{,}$(45)where the coefficients 0.48 and 0.5 are taken from Table 6 in Hollenbach & McKee (1989). The quantity G₀ represents the interstellar radiation field (IRF) normalized to the solar neighborhood value, i.e., we make here the usual assumption that the cosmic ray intensity scales with the IRF. We calculate its value as G₀ = SFR_loc/SFR_MW, where SFR_loc is the local SFR in our models and SFR_MW = 1.0 M_⊙ yr^-1 is the SFR in the Milky Way (Robitaille & Whitney 2010).

3.5.2. Photoelectric heating from small dust grains

The photoelectric emission from small grains induced by FUV photons can substantially contribute to the thermal balance of the interstellar medium at temperatures ≲ 10⁴ K. We use the results of Bakes & Tielens (1994) and express the photoelectric heating rate per unit volume as ${\mathrm{\Gamma }}_{\mathrm{ph}}\mathrm{=}\mathrm{1.0}\mathrm{\times }{\mathrm{10}}^{-24}{\mathit{ϵ}}_{\mathrm{ph}}\hspace{0.17em}{\mathit{G}}_{\mathrm{0}}{\mathit{n}}_{\mathrm{H}}{\mathit{\xi }}_{\mathrm{d}}\hspace{0.17em}\hspace{0.17em}\mathrm{[}{\mathrm{erg}\mathrm{cm}}^{-3}\hspace{0.17em}{\mathrm{s}}^{-1}\mathrm{]}\mathit{.}$(46)The photoelectric heating efficiency ϵ_ph is calculated from Fig. 13 of Bakes & Tielens. We also introduced the coefficient ξ_d = (M_d/M_HI)/(M_d/M_HI)_⊙ to take a (possibly) different dust to gas mass ratio (M_d/M_HI) in DGs into account. The latter is found using the following relation between the oxygen abundance (by number density) and the dust-to-gas mass ratio in dwarf galaxies (Lisenfeld & Ferrara 1998): $\mathrm{12}\mathrm{+}\log \left(\frac{\mathrm{O}}{\mathrm{H}}\right)\mathrm{=}\mathrm{0.52}\hspace{0.17em}{\mathit{C}}_{\mathrm{d}}\hspace{0.17em}\log \left(\frac{{\mathit{M}}_{\mathrm{d}}}{{\mathit{M}}_{\mathrm{HI}}}\right)\mathit{,}$(47)where the normalization constant C_d = −5.0 is derived from the corresponding quantities of the Small Magellanic Cloud. For the solar dust to gas mass ratio, we adopt (M_d/M_HI)_⊙ = 0.02, typical for the interstellar medium in the solar neighborhood.

4. Solution method

The equations of gas and stellar hydrodynamics are solved in cylindrical coordinates with the assumption of axial symmetry using a time-explicit (except for the gas internal energy update due to gas cooling/heating), operator-split approach similar in methodology to the ZEUS code (Stone & Norman 1992). We use the same finite-difference scheme but, contrary to the ZEUS code, we advect the internal energy density ϵ rather than the specific one ϵ/ρ_g because the latter may give rise to undesired oscillations in some test problems (e.g., Clarke 2010). Below, we briefly review the numerical scheme. The test results of our implementation of the code are provided in the Appendix.

The solution procedure consists of three main steps. Step 1: we update the gas and stellar components due to star formation and feedback. The corresponding ordinary differential equations describing the time evolution of ρ_g, ρ_i, ρ_gv_g, ϵ, ρ_s, ρ_sv_s, and Π_ij due to the source and sink terms involving , , , , and Γ_∗ in Eqs. (1)−(4) and (6)−(8) are solved using a first-order explicit finite-difference scheme. In Step 1, we also calculate the total gravitational potential of the gas and stellar components by solving for the Poisson Eq. (11) using the alternative direction implicit method described in Black & Bodenheimer (1975). The gravitational potential at the outer boundaries is calculated using a multipole expansion formula in spherical coordinates (Jackson 1975). To save on computational time, we invoke the gravitational potential solver only when the relative change in the total (gas plus stellar) density exceeds 10^-5 (see Stone & Norman 1992, for details).

Step 2: we solve for the gas hydrodynamics Eqs. (1)−(4), excluding the source and sink terms that have been taken into account in the previous step. The solution is split into the transport and source substeps. During the former, advection is calculated using a third-order accurate piecewise-parabolic scheme of Colella & Woodward (1984) modified for cylindrical coordinates as described in Blondin & Lufkin (1993). During the latter, we update the gas momenta ρ_gv_g due to gravitational and pressure forces and the gas internal energy ϵ due to compressional heating using a solution procedure described in Stone & Norman (1992).

At the end of Step 2, a fully implicit backward Euler scheme combined with Newton-Raphson iterations is used to advance ϵ in time due to volume cooling Λ and heating Γ rates. To monitor accuracy, the total change in ϵ in one time step is kept below 10%. If this condition is not met, we employ local subcycling in a particular cell by reducing the time step in the backward Euler scheme by a factor of two (as compared to the global hydrodynamics time step), and making two individual substeps in the cell where accuracy was violated. The local time step may be further divided by a factor of two until the desired accuracy is reached. The local subcycling help to greatly accelerate the numerical simulations as one need not to repeat the solution procedure for every grid cell.

In Step 2, we also add tensor artificial viscosity to Eqs. (3) and (4) to smooth out strong shocks as described in Stone & Norman (1992). As pointed out in Tasker et al. (2008) and discussed in Clarke (2010), the ZEUS code with the von Neumann & Richtmyer (1950) definition of the artificial viscosity performs poorly on multidimesional tests, such as the Sedov point explosion, yielding nonspherical solutions that may overshoot or undershoot the analytic solution. We demonstrate in Sect. A.2 that using the tensor artificial viscosity and the proper choice of the artificial viscosity parameter C₂ = 6 can greatly improve the code performance. However, contrary to the Stone & Norman suggestion to discard the off-diagonal terms of the viscous stress tensor, we demonstrate that in fact these terms should be retained and the most general formulation of the artificial viscosity stress tensor should be used, i.e., ${Q}\mathrm{=}\{\begin{array}{c}\\ {{l}}^{{2}}{{\rho }}_{{g}}{\nabla }{·}{{v}}_{{g}}\left[{\nabla }{{v}}_{{g}}{-}{e}{(}{\nabla }{·}{{v}}_{{g}}{)}{/}{3}\right]& \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{if} \mathrm{\nabla }\mathrm{·}{\mathit{v}}_{\mathrm{g}}\mathit{<}\mathrm{0}\\ \mathrm{0}& \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{if} \mathrm{\nabla }\mathrm{·}{\mathit{v}}_{\mathrm{g}}\mathrm{\ge }\mathrm{0}\mathit{,}\end{array}$(48)where e is the unit tensor, l = C₂max(dx), and dx is the grid resolution in each coordinate direction. For the reasons of the code stability, all components of the viscous stress tensor Q should be defined at the same position on the grid, i.e., at the zone centers (contrary to the ZEUS code, which defines the off-diagonal components of rank-two tensors at the zone corners). The Courant limitations on artificial viscosity may be substantial in regions with strong shocks. Therefore, we apply local subcycling (in the same manner as with the internal energy update due to heating/cooling) when the local viscous time step becomes a small fraction of the global hydrodynamics time step.

In Step 3, we solve for the stellar hydrodynamics Eqs. (6)−(8), excluding the source and sink terms that have been taken into account during Step 1. The solution procedure is similar to that for the gas hydrodynamics equations, since both systems are defined essentially on the same grid stencil, which is a powerful advantage of the stellar hydrodynamics approach. All components of the velocity dispersion tensor are defined at the cell centers to be consistent with the definition of the stellar artificial viscosity tensor Q. The latter is defined with an analogy to its gas counterpart with the gas density and velocity in Eq. (48) substituted with the stellar density and velocity. Additional terms −∇·Q and −Q_ik:(∇v_s)_kj − Q_jk:(∇v_s)_ki have been added to the r.h.s. of Eqs. (7) and (8), respectively, to take the viscous stress and heating due to artificial viscosity into account.

Finally, the global time step for the next loop of simulations is calculated using the Courant condition modified to take star formation and feedback into account. In particular, two time step delimiters of the form, $\begin{array}{ccc}\mathrm{△}{\mathit{t}}_{\mathrm{s}\mathit{.}\mathrm{f}\mathit{.}}& \mathit{<}& \frac{{\mathit{C}}_{\mathrm{3}}{\mathit{\rho }}_{\mathrm{g}}}{\mathrm{𝒟}}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{(}\mathrm{star}\mathrm{formation}\mathrm{)}\mathit{,}\\ \mathrm{△}{\mathit{t}}_{\mathrm{f}\mathit{.}\mathrm{b}\mathit{.}}& \mathit{<}& \frac{{\mathit{C}}_{\mathrm{3}}{\mathit{\rho }}_{\mathrm{s}}}{\mathrm{𝒮}}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{(}\mathrm{stellar}\mathrm{feedback}\mathrm{)}\end{array}$are added to the usual Courant condition to avoid obtaining negative values of the flow variables. The safety coefficient C₃ is set to 0.5. Some modification to the classic Courant condition is required for the stellar hydrodynamics part because the stellar velocity dispersion tensor Π is anisotropic in general. To calculate the stellar hydrodynamics time step, we use the following form: $\mathrm{△}{\mathit{t}}_{\mathrm{st}}\mathrm{=}\begin{array}{c}\min \end{array}\left[\frac{\mathrm{d}\mathit{x}}{\mathrm{\left|}{\mathit{v}}_{\mathrm{s}}\mathrm{\right|}\mathrm{+}\max {\overline{\mathrm{\Pi }}}_{\mathit{ii}}}\right]\mathit{,}$(51)where $\max {\overline{\mathrm{\Pi }}}_{\mathit{ii}}$ is the maximum diagonal element of the diagonalized velocity dispersion tensor and the minimum is taken over all grid zones.

5. The chemodynamical evolution of model dwarf galaxies

Our numerical code has been extensively tested using a standard set of test problems as described in the Appendix. However, even a perfect performance on the standardized tests does not guarantee meaningful results on real problems when all parts are invoked altogether.

Therefore, in this section, we consider the evolution of several representative models of DGs to test the overall code performance and to demonstrate the simulation methodology and recipes described in the previous sections. In particular, we wish to emphasize the effect of stellar dynamics and initial conditions on the simulation outcomes. Moreover, we show the results of simulations with different degrees of rotational support against gravity parameterized by the α-parameter defined as α = v_φ/v_circ, where v_φ and v_circ are the rotation and circular velocities, the latter being the maximum velocity calculated from the assumption that all support against gravity comes from rotation (i.e., zero gas pressure gradients). We consider three different degrees of rotational support: α = 0.25, 0.5, 0.9. The resulting initial geometries vary from flat-like for large values of α to roundish for small α.

Our initial configuration consists of a self-gravitating gaseous disk submerged into a fixed DM halo described by a cored isothermal sphere (e.g., Silich & Tenorio-Tagle 2001). We focus on models with a DM halo mass of 10⁹M_⊙ and initial gas temperatures of 10⁴ K. We apply the solution procedure described in Vorobyov et al. (2012) to construct self-gravitating equilibrium configurations for different values of α. The initial gas metallicity Z_g is set to 10^-2Z_⊙, the ratio of specific heats γ to 5/3, and the redshift is set to zero. The number of grid zones is 1000 in each coordinate direction, resulting in the effective numerical resolution of 8 pc everywhere throughout the computational domain. We impose the equatorial symmetry at the midplane and the reflection boundary condition at the axis of symmetry. At the other boundaries, a free outflow condition is applied.

Figure 3 presents the initial gas surface density Σ_g (left panel) and the Toomre Q-parameter (right panel) as a function of the midplane distance in the α = 0.25 model 1 (dash-dotted lines), α = 0.5 model 2 (dashed lines), and α = 0.9 model 3 (solid lines). Evidently, the α = 0.25 model is characterized by the most centrally condensed gas distribution. The horizontal dotted line marks a critical Q-parameter of 1.5 below which star formation is supposed to occur, according to the stability analysis of self-gravitating disks by Toomre (1964) and Polyachenko et al. (1997)⁵. Evidently, the α = 0.25 and α = 0.5 models are strongly inclined to star formation in the inner 0.5 kpc, while outside this radius the conditions are not favorable because of a rather steep drop in the gas surface density with distance. On the other hand, the α = 0.9 model is characterized by a much shallower profile of Σ_g, resulting in a significantly larger galactic volume prone to star formation. The parameters of our models are listed in Table 1, where M_g(r ≤ 1 kpc), M_g(r ≤ 3 kpc), and M_g(total) are the gas masses inside a spherical radius of 1 kpc, 3 kpc, and the total gas mass in the computational domain of 8.0 × 8.0 kpc².

Evidently, the gas mass increases with increasing α, which is explained by the fact that models with higher rotation can support a higher gas mass against gravity of the DM halo. The maximum gas mass is also limited by the action of gas self-gravity (Vorobyov et al. 2012). As the α = 0.5 model demonstrates, switching off the gas self-gravity leads to a factor of two increase in the total gas mass.

When constructing the initial gas density distribution in models 1, 2, and 3, we chose the same “seed” value for gas number density in the galactic center equal to 10 cm^-3. We did this to make the initial volume densities and temperatures in the galactic center similar in each model. We could have decreased the seed value in models 2 and 3, trying to adjust the total mass in these models to that obtained in model 1, but then the initial conditions for star formation would also have changed. In addition, it was not clear which mass to take as the reference value because different values are obtained when integrating over different radii because of variations in the radial density profiles.

5.1. The role of α

In Fig. 4 we show a comparison of the gas volume density (ρ_g) distribution in three model galaxies differing for the value of α. Four representative time snapshots are taken for each model. The α = 0.9 model 3 is characterized by the largest degree of flattening because of the large centrifugal forces. A notable flaring in the initial gas density distribution of model 3 can be seen in the upper right panel. This means that the vertical distribution at R = 0 is characterized by a very steep density and pressure gradient. The overpressurized superbubble created by the stellar feedback preferentially expands along this direction and a bipolar outflow is soon created. The external parts of the disk are not affected by the bipolar outflow and even after t = 150 Myr the gas density of the disk at radii larger than 1 kpc is quite large. Some of the material in the bipolar outflows falls back onto the disk. As a result, the fraction of pristine gas lost by this model is relatively small in spite of the fact that the galactic outflow is very prominent and affects a large fraction of the computational domain. However, the freshly produced metals can be easily channelled out of our model galaxy.

At the other extreme is the α = 0.25 model 1. Here, the initial gas density configuration is almost spherical and the pressure gradient is only slightly steeper in the vertical than in the horizontal direction. A bipolar outflow soon develops, but there is also significant gas transport along the disk and, after 100 Myr, all the gas in the central part of the galaxy (up to a radius of almost 5 kpc) has been completely blown away.

This effect is illustrates in Fig. 5, where the gas surface density as a function of radius is plotted for our three reference models at the same moments in time as in Fig. 4. After 100 Myr there is almost no gas left in the inner 3 kpc in the α = 0.25 model 1, whereas a central hole in the α = 0.9 model 3 is of a much smaller radius, hardly exceeding 0.2 kpc, and the gas distribution at R> 2 kpc is nearly undisturbed. This different galactic evolution, depending on the initial distribution of gas, has been already described in the literature (see in particular Recchi & Hensler 2013), but here we show it in a more realistic context of self-gravitating initial gas configurations.

A direct way to appreciate the role of α on the evolution of our model galaxies is to calculate the ejected gas mass for various models. Figure 6 presents the fraction of ejected gas mass vs. time in models 1–3. We find this value by calculating the amount of gas that leaves the computational domain (8 × 8 kpc²) with velocities greater than the escape velocity at the computational boundary. The fraction of ejected mass is very large in the α = 0.25 model because of the blow-away described above, but it is very small for the α = 0.9 model because only a small fraction of the gas originally present in the disk is involved in the galactic outflow. The dash double-dotted line presents the fraction of ejected mass in the α = 0.5 model 2, in which the motion of star was artificially suppressed (see Sect. 5.2 for details.)

The different gas dynamics of the three reference models shown in Fig. 4 also has dramatic consequences for the star formation histories. Figure 7 presents the SFR vs. time. Evidently, the strong star formation feedback in the α = 0.25 model shuts off star formation after ~30 Myr. The α = 0.5 model behaves similarly, although star formation is shut off slightly later. On the other hand, the α = 0.9 model is able to sustain star formation for a longer time, mainly because the disk still remains gas-rich and the conditions for the star formation to occur are fulfilled in many grid cells.

Figure 8 presents the distribution of the stellar volume density ρ_s as a function of time and for the same models and at the same time instances as in Fig. 4. In addition, the right-hand column presents the volume density of high-mass stars ${\mathit{\rho }}_{\mathrm{s}}^{\mathrm{h}}$. We can clearly see that the stellar distribution in the α = 0.25 model is more homogeneous and isotropic than in the α = 0.9 model, wherein many stars tend to be located along the polar axis. Owing to a preferable concentration of gas toward the midplane in the α = 0.9 model, however, a significant fraction of stars is also distributed near the midplane, constituting a stellar disk. This effect is much less evident for the models with α = 0.5 and α = 0.25. To quantify this, we calculated the fraction of the total stellar mass confined within a vertical height of 1 kpc and 2 kpc. It turns out that in the α = 0.25 model these values are 37% and 62%, respectively, whereas in the α = 0.9 model they are 90% and 94%, indicating a much more compact stellar distribution in the higher-α model. Finally, a comparison of the distribution between ρ_s (the total stellar density) and ${\mathit{\rho }}_{\mathrm{s}}^{\mathrm{h}}$ in the α = 0.9 model indicates that stars are preferentially forming at low galactic altitudes ≲ 1 kpc, but are later spread out owing to their own initial dispersion (5 km s^-1) and, especially, owing to significant bulk velocities inherited from parental gas clouds. The off-plane star formation is taking place predominantly in dense gaseous clumps and on the cavity walls.

5.2. The role of the stellar motion and the initial gas distribution

We employ the moments of the collisional Boltzmann equation to describe the motion of stars. In addition, we start our simulations from self-consistent initial equilibrium conditions that were built taking gas self-gravity into account. In some studies, the initial configuration is constructed without taking gas self-gravity into account, although during the dynamical evolution the gas self-gravity may be included. It is instructive to study how these two new recipes affect the dynamics of our model galaxies.

In Fig. 9 we show the gas volume density in three models with α = 0.5 at four representative evolution times. In particular, the left-hand column shows the evolution of our reference model 2, already shown in Fig. 4. The middle column shows model 4, in which the stellar dynamics has been artificially suppressed. Stars are allowed to form, but they are effectively motionless and are fixed to their birthplaces until they die. Evidently, the final gas distribution in model 4 is clumpy: the galactic volume is filled with moving cometary-like clouds. The increase in the number of clumps is caused by a reduced feedback in the absence of stellar motion. The moving clumps leave behind the newly born stars so that part of the stellar energy is released outside the clumps, which makes it easier for the clumps to survive. In the presence of stellar motion, the newly born stars follow the gas clumps (at least for some time) because they inherit gas velocity at the formation epoch, releasing the stellar energy inside the clumps and dispersing most of them.

The final gas distribution in model 2 is smooth and shows no cometary-like structures. At the same time, the outflow in model 4 clears a larger hole near the midplane than in the reference model. The combination of these two effects results in a rather chaotic gas surface density distribution shown in the middle panel of Fig. 5 with dotted lines. Model 4 is characterized by a somewhat higher star formation feedback, as is indicated by the fraction of ejected gas mass shown by the dot-dot-dashed line in Fig. 6.

The increased feedback in the model without stellar motions is likely because the energy of both supernova explosions and stellar winds is released in regions characterized by a systematically higher temperature and lower density than in the reference model. This is because star formation creates regions of slightly lower density compared to the surrounding gas, and stellar winds stir and heat up the circumstellar gas. If stars are motionless, then their feedback is confined to their birthplaces. On the other hand, if stellar motion is allowed, stars can move to neighboring regions, where the density is likely to be higher and where stellar winds have not yet heated the gas. The feedback will be thus less effective. This conclusion is at variance with that of Slyz (2007), who showed that stellar motion makes the feedback more effective and can even solve the overcooling problem. This different outcome is probably because of the lack of feedback from stellar winds in Slyz’s models.

Figure 10 presents the stellar volume density distributions for the same models as in Fig. 9. The left-hand column corresponds to the reference model 2, already shown in Fig. 8, while the middle column presents model 4, in which stellar motions are suppressed. Evidently, models with and without stellar motion produce strikingly different stellar distributions. The stellar distribution in the reference model is rather smooth, whereas in the model without stellar motions most of the stellar content is confined to the midplane and the inner region. In fact, the fraction of the total stellar mass confined within a vertical height of 0.3 kpc (from the midplane) is 82% for the reference model, but this value increases to 98% for the model without stellar motions. The curious stellar stripes visible in the middle column are caused by the development of a gaseous supershell: many stars tend to form on the walls of this supershell and, as a consequence, their position reflects the shell’s expansion. The lack of stellar motion makes these structures permanent. The fact that these structures are not observed in real galaxies outlines the importance of a proper treatment of stellar dynamics in numerical simulations of DGs.

Finally, we investigate the question of how the choice of the initial gas density distribution can influence the subsequent evolution of dwarf galaxies. For this purpose, we choose α = 0.5 and construct two initial gas density distributions: one taking gas self-gravity into account and the other without it. To compare the models with and without-self gravity easier, we require that the gas mass within the 3-kpc radius be similar in both models. The model with the initial equilibrium distribution taking gas self-gravity into account is in fact our model 2 considered above, while its counterpart is further referred as model 5 (see Table 1 for model parameters).

The initial radial profiles of the gas surface density and the Toomre Q-parameter in model 5 are shown in Fig. 3 by the dash-dot-dotted lines. Not surprisingly, the initial gas distribution in model 2 is more centrally condensed than in model 5: gas self-gravity creates an additional gravity pull toward the galactic center⁶, resulting in a denser initial configuration. For instance, the total gas mass in the inner sphere of 1 kpc in the initial model without self-gravity is ≈ 10⁷M_⊙, while this value rises by a factor of 2.8 in the model with self-gravity. As a consequence, the Toomre Q-parameter is initially greater than 2.0 everywhere and star formation is effectively suppressed.

The right-hand column in Fig. 9 presents the gas volume density distribution in model 5. We emphasize that the gas self-gravity in this model is only neglected when building the initial distribution; in the subsequent evolution the gas self-gravity is however turned on again. Evidently, the evolution of model 5 is notably slower owing to the fact that star formation is initially suppressed. The gas cools and contracts, and only after 120 Myr the conditions for star formation become favorable in the galactic center, whereas a vigourous star formation feedback is seen in model 2 already after 50 Myr of the evolution. The initial slow start in model 5 is consistent with a long dynamical timescale in this model. For the mean gas volume density in the inner 1 kpc of 4 × 10^-25 g cm^-3, the corresponding free-fall time in model 5 is approximately 110 Myr, which is consistent with the time epoch of the first star formation episode in this model. Because of significant thermal support and radially increasing dynamical timescale (density drops with radius), it takes several dynamical timescales to adjust to a state similar to that of model 2. By the end of simulations in model 5 (800 Myr), the outflow is still not as strong as in its counterpart model 2 and the central gap has not yet developed.

The right-hand column in Fig. 10 presents the stellar volume density distribution in model 5. A visual inspection of the l.h.s. and r.h.s. models reveals that the stellar distributions in models 2 and 5 are in general similar (but note the difference in ages 150 vs. 800 Myr), though the stars in the latter model are distributed somewhat more irregularly. By the end of numerical simulations, the total mass of formed stars in the 800-Myr-old model 5 is still a factor of 4 lower than in the 150 Myr-old model 2, highlighting a much slower evolution in the model that starts from a non-self-consistent initial configuration, i.e., without initial gas self-gravity.

6. Discussion and model caveats

We did not aim for a detailed comparison of our numerical simulations with observed DGs, which is our task for subsequent studies. Instead, we chose three typical models for DGs with different rotational support against gravity to test the code performance on global simulations, and not just on a predefined set of classic test problems. The current version of the code assumes axial symmetry and neglected chemical reactions, which leaves the door open for future improvements. At the same time, we managed to develop an efficient mathematical recipe for describing the phase transition from stars into gas using the continuity equation for the mean stellar age (see Sect. 2.6).

We chose some simulation parameters without a detailed parameter study or without an in-depth analysis. We avoid this kind of detailed study for the sake of clarity and to make the paper more concise and readable. A deeper parameter study will be performed in follow-up papers. It is nevertheless important to point out the following potential caveats in the choice of some parameters:

• ϵ_sw(see Sect. 2.7.3.): massive and intermediate-mass stars expel a significant fraction of their mass by stellar winds. Although their lifetimes shorten with mass, the mass-loss rates increase disproportionately, so that for solar abundances the total wind energy of most massive stars exceeds that of SNeII over the stellar lifetime. For radiation-driven HII regions, Lasker (1967) stated that the ISM is only heated a very low amount and derived analytically an efficiency of 1%, because the ionizing radiation energy is almost instantaneously lost by line emission and only the expansion of HII regions steers up the surrounding ISM. Although the additional wind power amounts to almost 10% of the radiation energy, it is expected to deposit a larger fraction as the thermal and kinetic ISM energies. Since in most chemodynamical simulations (Theis et al. 1992; Samland et al. 1997) the Ly continuum radiation is applied with an efficiency of 10^-3, ϵ_sw is assumed to 100%. This value is also used here.

• ϵ_SN(see Sect. 2.7.1): for the thermalization efficiency of the SNII explosion, i.e., the fraction of the explosion energy transferred to gas thermal energy, we have adopted the value ϵ_SN = 1.0 for simplicity and to avoid further parameter studies. The thermalization efficiency has been subject to considerable studies in the past. Thornton et al. (1998) performed single 1D supernova explosion models and derived ϵ_SN = 0.1. Estimates from cumulative SNII explosions, forming superbubbles, suggest a significant dependence of ϵ_SN on the temperature and density of the external medium (Recchi et al. 2001; Recchi 2014). A different approach, where the contribution of a whole stellar population is considered (Melioli & de Gouveia Dal Pino 2004) clearly shows that ϵ_SN is a function of time. We performed a test simulation adopting ϵ_SN = 0.1. The extent of the galactic winds was reduced in this case, but not by a large factor. In fact, a great deal of self-regulation is achieved with our formulation for the SFR (see Eq. (30)). Further, the larger ϵ_SN, the larger the heating of the surrounding medium, and the higher the probability that the star formation is quenched for a period of time (see Koeppen et al. 1995, for more details on self-regulated star formation). Thus, although we plan to make more careful selection of ϵ_SN in follow-up papers, this parameter does not appear change the results of the simulations drastically.

• n_cr (see Sect. 2.7.4): the threshold density above which star formation is allowed, n_cr = 1 cm^-3 has been chosen in compliance with other numerical studies (Springel & Hernquist 2003; Recchi et al. 2007). Also in this case, many different values for n_cr have been chosen in the literature and detailed studies on the effect of this parameters have been performed. In the framework of ΛCDM cosmological simulations, n_cr of 0.1 cm^-3 is generally used (Katz 1992) and its effect on the star formation was studied, e.g., by Kay et al. (2002), but also applied to much larger resolutions of smooth-particle hydrodynamics applications (Stinson et al. 2006). It appears that an extremely large value of (up to 100 cm^-3Governato et al. 2010) is required to reproduce the observed structural properties of DGs (but see also Pilkington et al. 2011). Quite generally, the star formation density threshold should be set based on the density at which gravitational instabilities can be resolved (see, e.g., Truelove et al. 1997). This leads, for instance, Stinson et al. (2013) to assume a threshold of 9.3 cm^-3 based on their adopted resolution. Clearly, our choice of n_cr is simplistic and deeper studies are required in follow-up publications. In our simulation, only a tiny fraction of gas reaches densities larger than a few tens of cm^-3, therefore the adoption of very large values of n_cr would lead to extremely localized and sporadic star formation. The very large feedback achieved in the simulations of Governato et al. (2010) is certainly able to produce modeled DGs with structural properties similar to those observed, but it is not clear whether the chemical properties can be equally well reproduced. A large amount of feedback can create galactic winds, which are able to expel the large majority of metals that are freshly produced in the galaxy, thus at the end of the evolution the model galaxy should be almost deprived of heavy elements.

• Λ_h (see Sect. 3.4): the cooling function typical of hot gas Λ_h, which has been used to calculate the cooling timescale of the hot SN ejecta, has been estimated based on a temperature of the hot ejecta of 10⁸ K. It must be admitted that this temperature only holds just after the SN explosion, but the ejecta cools significantly afterward. Nevertheless, the thermal SN energy has been varied artificially to probe its influence on galactic winds (Dalla Vecchia & Schaye 2012). However, we are in the cooling regime dominated by bremsstrahlung, and the temperature dependence of the cooling function is not so severe (the cooling is approximately proportional to $\sqrt{\mathit{T}}$ in this temperature range).

7. Conclusions

We presented the stellar hydrodynamics approach for modeling the evolution of DGs, paying special attention to the co-evolution of gas and stars. The distinctive feature of our models is how we treat the stellar component. Unlike many previous studies that use various versions of the N-body method to compute the dynamics of stars, we employ the stellar hydrodynamics approach based on the moments of the collisionless Boltzmann equation, originally introduced by Burkert & Hensler (1988), and applied and further developed by Theis et al. (1992), Samland et al. (1997), Vorobyov & Theis (2006, 2008), Mitchell et al. (2013), Kulikov (2014). For the first time, we provided an effective mathematical recipe for describing the death of stars in the framework of stellar hydrodynamics using the concept of a mean stellar age. We performed an extensive testing of our numerical scheme using classical and newly developed test problems, paying particular attention to the proper recipe for artificial viscosity to achieve the best agreement with the analytical solution for the Sedov test problem.

We demonstrate the success of our numerical approach by simulating the evolution of three DGs differing in the amount of rotation described by the spin parameter α = 0.25, α = 0.5, and α = 0.9. In agreement with previous studies, models with low values of α are characterized by stronger outflows, ejecting a larger amount of gas from the galaxy. This is because the initial gas distribution in the low-α models is almost spherical, so that the gas pushed by the stellar feedback expands almost isotropically. In contrast, models with high values of α exhibit a disk-like initial gas distribution with a pronounced radial flaring. The pressure gradient in the vertical direction in these models is very steep, which favors the development of bipolar galactic outflows. The transport of gas in the horizontal direction (i.e., along the galactic midplane) is very limited, and therefore this model is able to keep most of the initial gas bound to the galaxy. The fraction of ejected gas (i.e., the fraction of gas that leaves the computational domain with a speed higher than the escape speed) is only a few percent for the α = 0.9 model, but it raises to 80% for the model with α = 0.25.

At the same time, models with a low value of α produce a more diffuse stellar population than models with higher α. Only models with α ≥ 0.5 reveal the formation of a definite stellar disk component. We demonstrate the importance of the stellar dynamics by artificially turning off stellar motions. In this case, the resulting stellar population is mostly concentrated to the midplane and shows stripe-like structures emanating from the galactic center, which are not observed in real galaxies.

Stars start heating the natal gas almost instantaneously owing to the effect of stellar winds neglected in many recent hydrodynamical simulations. We found that if stellar dynamics is suppressed, Type II supernovae (the main source of mechanical feedback) explode in a medium heated and diluted by the stellar winds. In this kind of an environment, radiative losses are small and the supernova feedback is extremely effective. If, on the other hand, stellar motions are allowed, they can move from the natal site (where the gas has been heated and stirred by the effect of stellar winds) to neighboring regions before exploding as supernovae. These regions are characterized by lower temperatures and higher densities, so that the supernova feedback becomes less effective. It is thus crucial to correctly simulate stellar dynamics (together with all relevant sources of feedback) to study the interplay between stars and the interstellar medium.

Unlike many other studies of DGs, we used the self-consistent initial gas configuration that was constructed taking the gas self-gravity into account (Vorobyov et al. 2012). This gas configuration is usually more centrally condensed and prone to star formation than that constructed neglecting gas self-gravity. We compared the evolution of our model galaxy starting from these different initial configurations and found that models that start from non-self-gravitating initial configurations evolve much slower than those that start from self-gravitating configurations. In the former, a major episode of star formation and the development of prominent outflows may by delayed by 0.5−0.8 Gyr as compared to the latter. This is exclusively the effect of initial conditions (or, better to say, the effect of using non-self-consistent initial gas configurations), because the gas self-gravity is turned on in both models once the evolution has started.

In the future, a detailed comparison of our stellar hydrodynamics simulations with a more conventional N-body treatment of stellar dynamics is desirable to assess the weak and strong sides of the Boltzmann moment equation approach.

Online material

Appendix A: Testing gas hydrodynamics equations

Our numerical scheme for solving the equations of gas hydrodynamics (1)–(4) can be tested using a conventional set of test problems suitable for cylindrical coordinates. Below, we provide the results of six essential test problems, each designed to test the code performance at various physical circumstances. The code has also performed well on two additional tests: the two interacting blast waves and the Noh (or implosion) problems, but we do not provide the results for the sake of conciseness.

Appendix A.1: Sod shock-tube test

This test is often used to assess the ability of a numerical algorithm to accurately track the position of relatively weak shock waves and contact discontinuities. Initial conditions involve two discontinuous states along the z-axis, with a hot dense gas on the left and cold rarified gas on the right. More specifically, we set the pressure and density at z ∈ [ 0−0.5 ] to 1.0, while at z ∈ [ 0.5−1.0 ] the pressure is 0.1 and density is 0.125. The velocity of a γ = 1.4 gas is initially zero everywhere.

Filled circles in Fig. A.1 show the results of the Sod shock-tube test computed with a resolution of 200 grid zones and C₂ = 4 (the number of zones over which the shock is spread), while the solid lines plot the analytic solution at t = 0.235. It is evident that the code tracks well the positions of shock waves and contact discontinuities, which are resolved by 3–4 zones, in accordance with the chosen value of C₂. Decreasing C₂ results in sharper shocks and contact discontinuities but may increase the magnitude of anomalous spikes in the specific energy (bottom-right panel), and hence is not recommended. Increasing C₂ to 6 produces results similar to those for C₂ = 4, except that the shocks and contact discontinuities are now resolved by 5−6 zones.

Appendix A.2: Sedov point explosion

The Sedov explosion problem tests the code’s ability to deal with strong shocks. In contrast to the Sod shock-tube test, in which shock waves were of a pure planar geometry, this test involves a spherically symmetric strong shock wave. Hence, the point explosion is a powerful test on the code’s ability to render spherically symmetric structures on an essentially nonspherical grid stencil.

To initialize the explosion, we consider a cold (T = 10 K) homogeneous medium with number density n = 10^-2 cm^-3 and inject 10⁵¹ ergs of thermal energy (equivalent to one supernova) into the innermost grid cell. The adopted resolution is 300 × 300 grid cells and the size of each cell is 1.0 pc in each coordinate direction (z and r). We do not inject energy into a central sphere comprising a few cells near the coordinate origin, which is the usual practice to alleviate the problem of a nonspherical geometry. Instead, we consider the most difficult situation when energy is injected in just one, essentially cylindrical innermost cell.

Figure A.2 compares our derived density distributions along the z-axis (black line), the r-axis (pink line), and the diagonal z = r (blue line with circles) with the analytic solution given by the red line (Sedov 1959) at t = 1 Myr after the explosion. In particular, numerical simulations in the top left panel are obtained using the tensor artificial viscosity (AV) as defined in Sect. 4, while the density distribution in the top right panel is obtained using the “classic” scalar AV, as described by Richtmyer & Morton (1995). It is evident that the tensor AV yields a much better agreement with the analytic solution. The expanding shell has almost a perfect spherical shape (all three lines showing the numerical solution virtually coincide), notwithstanding the fact that the energy was injected into a cylindrically shaped innermost cell. On the other hand, the Richtmyer & Morton scalar AV yields a notable deviation from sphericity along the r and z axes, as seen in the top right panel of Fig. A.2. The Richtmyer & Morton solution also slightly overshoots the analytical solution. Regardless of the AV prescription used, the shock front is resolved by 2−3 grid zones.

We want to emphasize that all components of the artificial viscosity tensor Q should be used, including the off-diagonal terms. Neglecting the latter, as is done in, e.g., Stone & Norman (1992), leads to a significant deterioration of the numerical solution. The bottom right panel in Fig. A.2 demonstrates the effect of zeroing the off-diagonal terms in Q. As one can see, the numerical solution is skewed, notably undershooting the position of the shock front in the r-direction and along the diagonal direction but reproducing the shock front rather well in the z-direction.

A caution should be paid when choosing the AV softener C₂. As one can notice in Fig. A.2, the test was done with different values of C₂ for the tensor and scalar AV, 6 and 3, respectively. If we choose C₂ = 3 for the tensor AV, the position of the shock front notably undershoots the analytic solution, as is illustrated in the bottom left panel of Fig. A.2. Fortunately, once the proper value of C₂ is found, the tensor AV algorithm performs well regardless of the energy input, background density, or numerical resolution. We demonstrate this in Fig. A.3 showing the tensor AV performance for various choices of the input parameters as indicated in each panel and C₂ = 6. It is evident that the test results are virtually insensitive to a particular choice of initial conditions.

The sensitivity of the Sedov test to a particular choice of the AV softener C₂ is a well-known problem of numerical codes that evolve the internal energy (e.g., Tasker et al. 2008). Switching to the total energy helps to eliminate this problem (e.g., Clarke 2010). However, in numerical codes that make use of a staggered grid as our own, with vector and scalar variables defined at different positions on the grid, solving for the total energy equation is not a good option as it inevitably requires interpolating between the internal and kinetic energies to obtain the total energy (Clarke 2010) or applying corrections after each time step to conserve the total energy (Vshivkov et al. 2011). In either case, this practice, according to our experience, often comes at a cost of degraded performance on other test problems and is not recommended. Instead, we show that a proper choice of the artificial viscosity softener can help to resolve this problem without resorting to the total energy equation.

Appendix A.3: Collapse of a pressure-free sphere

The gravitational collapse of a pressure-free sphere is used to assess the code’s ability to accurately treat converging spherical flows in cylindrical coordinates. This test is also useful for estimating the performance of the gravitational potential solver on dynamical problems (our solver does well on the static configurations such as spheres and ellipsoids). To run this test, we set a cold homogeneous sphere of unit radius and density (for convenience, the gravitational constant is also set to unity) and let it collapse under its own gravity. The analytic solution to this problem only exists in the limit of an infinite cloud radius, describing the collapse of every mass shell (Hunter 1962). In the cylindrical geometry, however, we consider a cloud of finite radius with a sharp outer boundary to preserve the cloud sphericity. As a result, a rarefaction wave develops after the onset of the collapse, propagating toward the coordinate origin and necessitating complicated corrections to the analytic solution (Truelove et al. 1998).

The left panel in Fig. A.4 compares the results of our numerical simulations with the “uncorrected” analytic solution of (Hunter 1962, red line) at t = 0.535 (or 0.985 that of the free-fall time) when the initial density has increased by nearly three orders of magnitude. As in the previous test, we plot the numerically derived density along the z-axis (black line), r-axis (green line), and the z = r diagonal (blue line with circles). The numerical resolution is 200 × 200 grid zones. The cloud’s outer boundary is somewhat smeared out because of the action of the rarefaction wave. The peak density is also about 1% smaller than that predicted from the analytic solution. On the other hand, the code performs well on preserving sphericity of the collapsing cloud. In the right panel, we turn off the volume averaging corrections applied to the advection of the r-momentum in the r-direction as described in Blondin & Lufkin (1993), which results in the development of an artificial density spike near the z-axis and also in the distortion of a spherically symmetric collapse. This test demonstrates the utility of the volume averaging corrections introduced originally by Blondin & Lufkin (1993) in the PPA advection scheme in non-Cartesian coordinates.

Appendix A.4: Collapse of a rotating sphere

This test is invoked to assess the code’s ability to conserve angular momentum. The setup consists of a isothermal (T = 10 K) homogeneous sphere of unit radius and density, which rotates at a constant angular velocity Ω₀. The latter is found from the requirement that the ratio β of the rotational energy ${\mathit{E}}_{\mathrm{rot}}\mathrm{=}\mathit{I}{\mathrm{\Omega }}_{\mathrm{0}}^{\mathrm{2}}\mathit{/}\mathrm{2}$ to gravitational energy ${\mathit{E}}_{\mathrm{grav}}\mathrm{=}\mathrm{3}\mathit{G}{\mathit{M}}_{\mathrm{0}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}\mathrm{5}\mathit{R}\mathrm{\right)}$ is equal to 5%. Here, I = 2M₀R²/ 5 is the moment of inertia of a homogeneous sphere of radius R and mass M₀. The adopted value of β = 0.05 is close to an upper limit found in rotating prestellar molecular cloud cores and DM halos.

In the absence of any mechanisms for angular momentum redistribution, the mass with specific angular momentum less than or equal to K = ωr², where ω is the angular velocity at a distance r from the z-axis, should be conserved and equal to (Norman et al. 1980), $\mathit{M}\mathrm{\left(}\mathit{K}\mathrm{\right)}\mathrm{=}{\mathit{M}}_{\mathrm{0}}\left[\mathrm{1}\mathrm{-}{\left(\mathrm{1}\mathrm{-}\frac{\mathit{K}}{{\mathrm{\Omega }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{2}}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}\right]\mathit{.}$(A.1)A deviation from Eq. (A.1) would manifest the nonconservation of angular momentum in the numerical algorithm.

The solid line in Fig. A.5 shows the initial mass spectrum M(K), whereas the open circles present the mass spectrum obtained at time 1.03 × t_ff = 0.551, where t_ff is the free-fall time for a nonrotating sphere of the same density and radius. The numerical resolution is 200 × 200 grid zones. By t = 0.551, the density near the center of the sphere has increased by more than three orders of magnitude, but the deviation of the mass spectrum from the initial configuration is negligible everywhere except at lowest values of K. This deviation is a manifestation of angular momentum nonconservation near the rotational axis caused by imperfections in the interpolation procedure across the inner r-boundary. Nevertheless, the total mass of gas that suffers form the angular momentum nonconservation is just below 1%, and therefore should not affect notably the global evolution.

Appendix A.5: Kelvin-Helmholtz instability

The Kelvin-Helmholtz instability (KHI) occurs at the interface between two fluid moving with different velocities. In the context of this paper, this instability operates at the interface between the hot supernova ejecta and the cold swept-up shell and may cause an accelerated disintegration of the latter (e.g., Vorobyov & Basu 2005). Therefore, it is essential that the code be able to demonstrate the development of the KHI as expected from the linear perturbation theory.

The initial setup involves two fluids moving in opposite directions along the z-axis with the relative velocity v_z,rel = 0.2c_s, where c_s = 0.06 is the sound speed. The interface between the fluids is located at r = 0.5 and the densities of the inner and outer fluids are ρ_in = 1.0 and ρ_out = 0.1, respectively. We use periodic boundary conditions and the resolution is 200 × 200 grid zones on the 1.0 × 1.0 computational domain.

To trigger the instability, we impose a sinusoidal perturbation on the radial velocity v_r of the form ${\mathit{v}}_{\mathrm{r}}\mathrm{=}\mathit{\delta }{\mathit{v}}_{\mathrm{r}}\sin \mathrm{(}\mathrm{2}\mathit{\pi z}\mathit{/}\mathit{\lambda }\mathrm{)}\mathit{,}$(A.2)where the amplitude and wavelength of the perturbation are δv_r = v_z,rel/ 100 and λ = 1/6, respectively. The perturbation is applied to five neighboring grid zones on both sides of the fluid interface.

Figure A.6 demonstrates the growth of the KHI with time. There are six growing vortices, in agreement with the wavelength of the perturbation λ = 1/6, and the growth time is approximately 50−60, again in agreement with the characteristic growth time of the KHI for our choice of the fluid density and velocity, ${\mathit{\tau }}_{\mathrm{KH}}\mathrm{=}\frac{\mathit{\lambda }\mathrm{(}{\mathit{\rho }}_{\mathrm{in}}\mathrm{+}{\mathit{\rho }}_{\mathrm{out}}\mathrm{)}}{{\mathit{v}}_{\mathrm{z}\mathit{,}\mathrm{rel}}\sqrt{{\mathit{\rho }}_{\mathrm{in}}{\mathit{\rho }}_{\mathrm{out}}}}\mathrm{=}\mathrm{48.3.}$(A.3)We should note here that, for this particular test, gravity can be neglected. Therefore, all modes of perturbation, irrespective of their wavelengths (or wavenumbers) are unstable (see Vieser & Hensler 2007, for stability criteria when gravity is taken into account).

Appendix A.6: Overstability of a radiative shock

A useful test problem for hydrodynamic codes with optically thin radiative cooling is described in Stone & Norman (1993). The test involves setting a unform flow of gas with velocity v₀ against a reflecting wall. An adiabatic reflection shocks forms immediately at the wall and propagates upstream with velocity v_sh = v₀/ 3. After approximately one cooling time, the shock loses its pressure support and is then advected back to the wall at v_sh ≲ −v₀. After reflecting off the wall, the shock again becomes adiabatic and propagates outward to repeat the cycle.

Figure A.7 presents the test results for the cooling function with solar metallicity and ionization fraction f_ion = 10^-3. Gas heating is turned off for this test. The initial setup is identical to that described in Stone & Norman (1993), i.e., the initial gas number density and temperature are set to 10 cm^-3 and 10⁴ K, respectively. The position of the shock front for three flow velocities of 120 km s^-1, 130 km s^-1 and 150 km s^-1 are plotted with the solid, dash-dotted, and dashed lines, respectively. The shock position demonstrates regular pulses, as expected. The amplitude and period of the observed oscillations are somewhat different from those presented in Stone & Norman (1993), which simply reflect the difference in the adopted cooling functions. The upstream velocity of the adiabatic shock is close to the expected value, v₀/ 3. For instance, in the case of v₀ = 150 km s^-1, the maximum upstream velocity of the shock is 47 km s^-1. This test demonstrates the reliability of our solution scheme for the gas cooling and heating described in Sect. 3.

Appendix B: Testing stellar hydrodynamics equations

Testing the stellar hydrodynamics Eqs. (6)−(8) presents a certain challenge since stellar systems are described by a stellar dispersion tensor, which may be anisotropic, rather than by an isotropic pressure. Fortunately, some of the test problems considered above can be adapted to stellar hydrodynamics as well. The matter is that the equations of stellar hydrodynamics can be formally made identical to those of gas hydrodynamics for a specific set of initial conditions (note that the continuity equations are always identical). Indeed, for a one-dimensional flow along the z-direction (with zero gravity and no star formation), Eqs. (7) and (8) become $\begin{array}{ccc}& & \frac{\mathit{\partial }}{\mathit{\partial t}}\left({\mathit{\rho }}_{\mathrm{s}}{\mathit{v}}_{\mathrm{s}\mathit{,z}}\right)\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial z}}\left({\mathit{\rho }}_{\mathrm{s}}{\mathit{v}}_{\mathrm{s}\mathit{,z}}\mathrm{·}{\mathit{v}}_{\mathrm{s}\mathit{,z}}\right)\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial z}}\mathrm{(}{\mathit{\rho }}_{\mathrm{s}}{{\mathit{\sigma }}_{\mathit{zz}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{(}{\mathit{\rho }}_{\mathrm{s}}{{\mathit{\sigma }}_{\mathit{zz}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial z}}\mathrm{(}{\mathit{\rho }}_{\mathrm{s}}{\mathit{\sigma }}_{\mathit{zz}}^{\mathrm{2}}\mathrm{·}{\mathit{v}}_{\mathrm{s}\mathit{,z}}\mathrm{)}\mathrm{+}\mathrm{2}{\mathit{\rho }}_{\mathrm{s}}{\mathit{\sigma }}_{\mathit{zz}}^{\mathrm{2}}\frac{\mathit{\partial }{\mathit{v}}_{\mathrm{s}\mathit{,z}}}{\mathit{\partial z}}\mathrm{=}\mathrm{0.}\end{array}$These equations become identical to the corresponding gas dynamics equations for ρ_gv_z and ϵ, if we set ${\mathit{\rho }}_{\mathrm{s}}{\mathit{\sigma }}_{\mathit{zz}}^{\mathrm{2}}\mathrm{=}\mathit{P}$, $\mathit{ϵ}\mathrm{=}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}{\mathit{\rho }}_{\mathrm{s}}{\mathit{\sigma }}_{\mathit{zz}}^{\mathrm{2}}$, and γ = 3 (see Mitchell et al. 2013, for details). This fact allows us to directly compare the performance of both stellar and gas hydrodynamics code on some test problems considered in Sect. A or use analytic solutions available from the gas dynamics.

Appendix B.1: Sod shock-tube test

The initial setup for this test is identical to that considered in Sect. A.1. We run the test along the z-axis and set ${\mathit{\rho }}_{\mathrm{s}}{\mathit{\sigma }}_{\mathit{zz}}^{\mathrm{2}}$ and ρ_s at z ∈ [ 0−0.5 ] to 1.0, while at z ∈ [ 0.5−1.0 ] the z-component of the stress tensor is 0.1 and stellar density is 0.125.

Our numerical results for C_2s = 4 and t = 0.14 (open circles) are compared against the analytical solution (solid lines) in Fig. B.1. It is seen that our numerical scheme correctly reproduces the position of shock waves and contact discontinuities in the stellar fluid. On the other hand, small-scale oscillations are evident at shock positions, perhaps reflecting the lack of energy dissipation due to the collisionless nature of stellar hydrodynamics. An increase in the coefficient of artificial viscosity helps to reduce the amplitude of spurious oscillations but simultaneously increases the shock smearing and hence is only recommended if these oscillations lead to numerical instabilities.

Appendix B.2: Collapse of a cold sphere and angular momentum conservation

The collapse of a cold stellar sphere, both with and without rotation, produces very similar results to those presented in Sects. A.3 and A.4 and will not be repeated here. This ensures that our stellar hydrodynamics part is expected to perform well on problems involving gravitational contraction/expansion with and without rotation.

Appendix C: Testing the chemical evolution routines

It is currently common to find hydrodynamical codes coupled with passive scalars (or other numerical techniques), which are able to follow the evolution of the metals in a simulation. However, usually no test is performed to check the numerical accuracy of these routines. Actually, since a very vast literature on analytical solutions of the chemical evolution of galaxies and other astronomical objects exists (see, e.g., Tinsley 1980; Matteucci 2001; Hensler & Recchi 2010, and references therein) it is indeed not difficult to benchmark the chemical evolution routines of a chemodynamical code.

The “closed box” model of chemical evolution is based on the following assumptions: (i) the system is uniform and closed (there are no inflows or outflows); (ii) the IMF does not change with time; (iii) the gas is well mixed at any time; and (iv) stars more massive than a certain threshold mass M_thr die instantaneously whereas stars less massive than M_thr live forever. This last assumption, although very restrictive, is fulfilled if we look at specific chemical elements, in particular O and the other α-elements, since they are mostly synthesized by Type II SNe and the (long-living) low- and intermediate-mass stars contribute negligibly to their production.

In a system made up of only gas and stars (DM does not play any role here) we define the quantities, $\begin{array}{ccc}& & \mathit{R}\mathrm{=}{\mathrm{\int }}_{{\mathit{M}}_{\mathrm{thr}}}^{\mathrm{\infty }}\left[\mathit{M}\mathrm{-}{\mathit{M}}_{\mathrm{rem}}\mathrm{\left(}\mathit{M}\mathrm{\right)}\right]\mathit{\phi }\mathrm{\left(}\mathit{M}\mathrm{\right)}\mathrm{d}\mathit{M,}\\ & & {\mathit{y}}_{\mathrm{O}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{-}\mathit{R}}{\mathrm{\int }}_{{\mathit{M}}_{\mathrm{thr}}}^{\mathrm{\infty }}\mathit{M}{\mathit{p}}_{\mathrm{O}}\mathrm{\left(}\mathit{M}\mathrm{\right)}\mathit{\phi }\mathrm{\left(}\mathit{M}\mathrm{\right)}\mathrm{d}\mathit{M,}\end{array}$where M_tot = M_gas + M_stars, M_rem(M) is the remnant mass after a star of mass M has died and Mp_O(M) is the oxygen mass newly synthesized by the star of mass M (p_O(M) is the true yield of O). Under the initial condition M_gasf(0) = Const. = M_tot, M_stars(0) = 0, X_O(0) = X₀ we find that ${\mathit{X}}_{\mathrm{O}}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{O}}}{{\mathit{M}}_{\mathrm{gas}}}\mathrm{=}{\mathit{X}}_{\mathrm{0}}\mathrm{+}{\mathit{y}}_{\mathrm{O}}\ln \mathrm{\left(}{\mathit{\mu }}^{-1}\mathrm{\right)}\mathit{,}$(C.3)

where μ = M_gas/M_tot. This is the famous analytical solution of the closed box model, which only depends on the gas fraction μ.

We have integrated numerically Eqs. (C.1) and (C.2) by using the values of M_rem(M) and Mp_O(M) tabulated by Woosley & Weaver (1995), using Z = 0.01Z_⊙ as initial metallicity. We do not use Z = 0 as initial metallicity for two reasons: firstly, for Z(0) = 0 the cooling is very low and the resulting star formation is too weak; secondly, the yields p_O(M) tabulated by Woosley & Weaver (1995) at this metallicity oscillate considerably and that makes the calculation of y_O by means of Eq. (C.2) less accurate. By properly choosing a reference volume in the computational box (namely a volume encompassing all the metals produced and advected within a time span of 100 Myr), we can thus calculate the variation of μ with time and, from it, the analytical evolution of the oxygen mass fraction X_O with time. This is plotted in Fig. C.1 (solid line), together with the direct numerical calculation of M_O/M_gas, averaged over the same volume (dashed line). The two curves only converge after ~30 Myr. This is expected because we have to ensure that Type II SNe of all initial masses contribute to the production of O and ~30 Myr is approximately the lifetime of the less massive SNeII. After a time of ~40 Myr the two curves almost overlap, demonstrating that our chemical routines are able to accurately follow the overall chemical evolution of a galaxy.

---

Acknowledgments

This project was partly supported by the Russian Ministry of Education and Science Grant 3.961.2014/K and a Lise Meitner Fellowship, project number M 1255-N16 and RFBR Grant 15-01-00508. The simulations were performed on the Vienna Scientific Cluster (VSC-2). This publication is supported by the Austrian Science Fund (FWF). We are thankful to the anonymous referee for a very insightful review that helped improve the manuscript.

References

All Tables

All Figures