Osaka Feedback Model II: Modeling Supernova Feedback Based on High-Resolution Simulations

When multiple SN explosions occur clustered in space and time in stellar clusters, superbubbles are formed with the shell-formation time being a critical timescale for their growth. In Section 2.1.1, we analytically calculate the shell-formation time and investigate its dependence on metallicity.

We carried out 3D hydrodynamical simulations with the Athena++ code²²2https://www.athena-astro.app (Stone et al., 2020) with the second-order accurate van Leer time integrator, the HLLC solver, and second-order spatial reconstruction. We solved for the equations of continuity, momentum conservation, and energy while incorporating cooling, heating, and thermal conduction:

where $E$ is the energy density. We did not consider self-gravity and magnetic fields in our simulation. Assuming a mean molecular weight of $\mu=1.4$, the gas temperature is taken at $T=P/(1.1n_{\text{H}}k_{\rm B})$, where the hydrogen number density is $n_{\text{H}}=\rho/(1.4m_{\rm H})$. We followed Kim & Ostriker (2015) for the implementation of radiative cooling and heating. For the cooling function $\Lambda(T)$ at low ($T<10^{4}$ K) and high ($T>10^{4}$ K) temperatures of gas, we adopted the fitting formulae from Koyama & Inutsuka (2002)³³3The cooling function given by Eq. (4) of Koyama & Inutsuka (2002) contains two typographical errors. See Eq. (4) of Vázquez-Semadeni et al. (2007) for the corrected functional form. and piecewise power-law fits to the cooling function of Sutherland & Dopita (1993) with collisional ionization equilibrium⁴⁴4https://www.mso.anu.edu.au/~ralph/data/cool/, respectively, as shown in Fig. 1. Heating was applied, with the following heating function, only at $T<10^{4}$ K to model the photoelectric heating of warm/cold ISM:

We used different cooling functions at $T>10^{4}$ K for different metallicities, but we did not take into account the metallicity dependence at $T<10^{4}$ K. We did not consider the metallicity dependence of the heating function. At $Z\ll 1$, the low-temperature cooling and heating rates would be quite different due to reduced photoelectric heating and cooling by fine-structure C and O lines when the metal abundance is low. However, these differences do not affect the results of our simulations, because the shell formation time is determined by the cooling time of shock-heated gas at high temperatures. We imposed a temperature floor of $10^{3}$ K.

We followed El-Badry et al. (2019) to implement thermal conduction. For thermal conductivity $\kappa$ at low ($T\leq 6.6\times 10^{4}{\rm K}$) and high $(T>6.6\times 10^{4}{\rm K})$ temperatures of gas, we applied, respectively, the thermal conductivities due to the neutral atomic collision (Parker, 1963):

where $T_{4}=T/10^{4}\,{\rm K}$, and hydrogen plasma (Spitzer, 1962):

where $T_{7}=T/10^{7}\,{\rm K},\ n_{e,-2}=n_{e}/10^{-2}\rm cm^{-3}$ and $n_{e}$ denotes the electron number density. The conductive heat flux $q=\kappa\nabla T$ is limited by the energy flux that can actually be transported by the electrons $q_{\rm max}=\frac{3}{2}\rho c_{\rm s,iso}^{3}$ (Parker, 1963; Cowie & McKee, 1977). The effective thermal conductivity $\kappa_{\rm eff}$, which includes the saturation resulting from limiting the heat flux, becomes

To save calculation time, we imposed a ceiling on the thermal conductivity at the following value:

This value is sufficiently high and has negligible effects on our results. We note that our resolution is not high enough to resolve the conductive interface.

We ran 36 simulations by varying the initial values of the parameters as listed in Table 1. The spatial resolution was fixed in time and uniform over the simulation box with periodic boundary conditions (we did not use static mesh refinement (SMR) or adaptive mesh refinement (AMR)). The spatial resolutions were $\Delta x$ = 6, 3, and 0.75 pc for number densities of $n_{\text{H}}=$ 0.1, 1, and 10 ${\rm cm^{-3}}$, respectively, to satisfy the convergence condition by Kim & Ostriker (2015), $\Delta x<R_{\rm sf}/10$. To set up a turbulent initial condition, we start from a uniform density field with pressure of $2.0\times 10^{3}\,(n_{\text{H}}/1\,{\rm cm^{-3}})\,{k_{\rm B}}\,{\rm K}\,{\rm cm^{-3}}$. We then generate the decaying turbulence with only initial forcing and evolve to 1 Myr. The initial forcing field was normalized to a Kolmogorov (1941) power spectrum $E(k)\propto k^{-5/3}$, while the range of driving was set to $2\leq kL/2\pi\leq 20$. The velocity dispersion in the box was set to $\sigma=5$ km s^-1 adopted based on Larson’s law (Larson, 1981; Ohlin et al., 2019). We note that this velocity dispersion of 5 ${\rm km\,s^{-1}}$ is lower than observed values at the scale of the simulation box, and that the initial evolution of 1 Myr is not long enough for multiphase ISM to develop via thermal instability. As a result, the inhomogeneity in the background is very weak, and the bubble expansion can be much more spherical than in realistic cases, which in turn affects the development of the Kelvin–Helmholtz instabilities that are normally driven by the shear at interfaces.

Thermal energy of $10^{51}$ erg per SN was injected into cells whose centers were at a distance $<r_{\rm init}$ from the site of the SN explosion, where $r_{\rm init}$ is the radius within which the mass is 1 ${\rm M_{\odot}}$. When thermal energy was injected, a mass of 1 ${\rm M_{\odot}}$ was also simultaneously injected. We repeated this ten times, at intervals of $\Delta t_{\rm SN}$, setting $t=0$ at the time of the initial energy injection.

To study the time evolution of the physical quantities of the superbubble, we define the bubble radius $R_{\rm bub}$ as the largest radius which satisfies following two criteria: (i) $v_{r}>3\,{\rm km\,s^{-1}}$, and (ii) $t_{\rm dens}<4t_{\rm sc}$, where $v_{r}$ is the radial velocity averaged over a spherical shell of radius $r$, $t_{\rm dens}=\rho/\dot{\rho}$ indicates the timescale of the density change, and $t_{\rm sc}=\Delta x/c_{s}$ depicts the sound-crossing time. In practice, we compute $v_{r}$, $t_{\rm dens}$, and $t_{\rm sc}$ for each cell, draw averaged radial profiles of these quantities, and then apply the above criteria. We have described the time evolution of total mass, energy, and momentum within $R_{\rm bub}$ in Section 2.3.1.

To discuss the time evolution of the hot bubble, we also define the hot-gas radius as $R_{\rm hot}=(V_{\rm hot}/(4/3)\pi)^{1/3}$, where $V_{\rm hot}$ is the region at which $T>10^{4}$ K. The time evolution of $R_{\rm hot}$ is discussed in Section 2.4.

In this section, we derive the scaling relations for the time evolution of superbubble momentum and radius for the application to the SN feedback model in galaxy simulations.

To discuss the environmental dependence of the superbubble in a unified manner, we normalize each physical quantity by its value at the shell-formation time. The same was done for single SN explosion simulations by Kim & Ostriker (2015), and also by Kim et al. (2017) for only the time variable for multiple SN explosion simulations.

Figure 8 shows the time evolution of the radius, normalized by its value at the shell-formation time. The time evolution of the radius shows some deviations, but in all cases, we can see that the radius evolves according to

as was also found by Kim et al. (2017). When $t<t_{\rm sf,m}$, the superbubble is in an adiabatic period and expected to evolve with $R\propto t^{3/5}$ (Eq. 11). The good fit of the straight line with a slope of 0.5 can be attributed to the fact that the energy injection is discrete rather than continuous, resulting in an intermediate time evolution with the Sedov–Taylor solution $R\propto t^{2/5}$ (Eq. A4). The asymptote of the fitting line is due to the transition from the Sedov–Taylor solution to the adiabatic solution of the superbubble. Time evolution after shell-formation can also be approximated by a straight line with a slope of 0.5. The asymptotic behavior of the blue lines from the bottom is because the hot-gas radius $R_{\rm hot}$ could not correctly follow the radius of the inside of the shell when it was cooled down to $T<2\times 10^{4}\,{\rm K}$. The scaling relation $R\propto t^{1/2}$ agrees with the prediction from the momentum-driven model and not with that from the pressure-driven model (Appendix B). Substituting the values of $R_{\rm sf,m},\,t_{\rm sf,m}$ (Eqs. (13), (12)), the time evolution of radius and velocity are obtained as

Next, we show the time evolution of momentum in Figure 9. It shows two different trends after the shell-formation time. If the interval between SN explosions is longer than the duration of the pressure-driven snowplow (PDS) phase of a single SN explosion, $t_{\rm PDS}$ (Eq. A15), energy injection is considered to be discrete and the superbubble is expected to display an evolution over time that is different from the continuous case. Since thermal transfer to the shell was not considered in the derivation of $t_{\rm PDS}$, the actual duration of the PDS phase will be shorter than $t_{\rm PDS}$.

Here we find that the following fitting functions can describe the two separate regimes continuously as shown in Fig. 9:

where $\tau=t/t_{\rm sf,m}$. Since the momentum can be estimated as $p=(4/3)\pi R^{3}\rho\dot{R}$, when the time evolution of the radius is expressed as $R\propto t^{\alpha}$, $p\propto t^{4\alpha-1}$. Before shell formation, the superbubble grows adiabatically, and the time evolution of its radius is $R\propto t^{3/5}$; thus, the time evolution of momentum is $p\propto t^{7/5}$.

When $t\gg t_{\rm sf,m}$, there occurs a deviation from the fitting line. This was also seen when $\Delta t_{\rm SN}=1\,{\rm Myr},\,n_{\text{H}}=10\,{\rm cm^{-3}}$ (Fig. 7), because the shell velocity dropped below $3\,{\rm km\,s^{-1}}$ and its time evolution could not be followed.

In passing, we note that El-Badry et al. (2019) derived a criterion to distinguish discrete energy injection versus continuous limit in their Eq. (53). They derived this criteria only by considering the blast waves from individual SNe before reaching the shell, and their simulation results asymptote to their modified energy-driven solution with cooling at $t>t_{\rm subsonic}$. We only simulated ten SN explosions, while El-Badry et al. (2019) examined the evolution all the way to the continuous energy injection limit (although with 1-D simulation); therefore our time-scale shown in Fig. 9 is somewhat of shorter time range than that of El-Badry et al. (2019)’s work (see their Fig. 6). The subsonic timescale $t_{\rm subsonic}$ in Eq.(53) of El-Badry et al. (2019) is much longer than our $t_{\rm PDS}$, and therefore it is not straightforward to compare our results using their criteria based on $t_{\rm subsonic}$.

From the fitting function in Sec. 2.4, we derive the expression for momentum of a superbubble per SN, averaged over an initial stellar cluster mass function (ICMF). We assume the ICMF to be described by the power-law $dN/dM_{\mathrm{c}}\propto M_{\mathrm{c}}^{-2}$ with high-mass cutoff, where $M_{\mathrm{c}}$ is the mass of the stellar cluster (Krumholz et al., 2019) and stellar initial mass function (IMF) is fixed. We also assume SNe occurs at a constant rate. The ICMF can be translated to an SN number function (probability distribution function of the number of SN explosions occurring in a stellar cluster)

where $N_{\mathrm{SN}}$ is the number of SN explosions occurring in a stellar cluster and $A$ is the normalization factor. The range of the number of SN explosions is assumed as $N_{\mathrm{SN}}=N_{\rm SN,min}-N_{\rm SN,max}$, where we adopt the range of ($N_{\rm SN,min}$, $N_{\rm SN,max}$) = (1, 500), and the normalization factor is $A=1/\ln{500}$. The value of $N_{\mathrm{SN,max}}=500$ is the number of SN explosions expected to happen in the largest young massive star clusters (YMCs) in the Milky way with a mass of $M=5\times 10^{4}\,{\rm M_{\odot}}$ (Portegies Zwart et al., 2010). The choice of high-mass cutoff has little effect on the result since such high-mass YMCs are rare objects. The SNe interval is obtained by dividing their duration by their number, $\Delta t_{\mathrm{SN}}=(t_{\mathrm{SN,end}}-t_{\mathrm{SN,start}})/N_{\mathrm{SN}}$, where $t_{\mathrm{SN,start}},\,t_{\mathrm{SN,end}}$ are the times at the first and the last Type II SN after the formation of the stellar cluster, i.e., the minimum and the maximum lifetimes of stars causes Type II SNe. The lifetimes of stars depend on their mass and metallicity. We calculate SNe duration using CELib⁵⁵5https://bitbucket.org/tsaitoh/celib (Saitoh, 2017), which compute it using the metallicity-dependent stellar lifetime table by Portinari et al. (1998). We assume the mass range of the progenitors of Type II SNe to be 13–40 ${\rm M_{\odot}}$. In the range of $Z=10^{-6}$–$10^{-2}$, SNe duration is

Single SN explosion energy is set to $10^{51}$ erg. The averaged momentum per SN is

The final integral is well-fitted by the following function:

and the momentum input by ${\cal N}_{\rm SN}$ SNe (superbubble momentum, $p_{\rm SB}$) is estimated as

Here we used a slightly different font of ${\cal N}_{\rm SN}$ for the total number of SNe for a star particle in the simulation (i.e., a collection of stellar clusters, rather than a population of stars in a single stellar cluster).

We similarly calculate the shock radius of SN feedback of a star particle. Assuming that stellar clusters follow the power-law ICMF, we determine the shock radius of the SN feedback as the radius of a sphere whose volume is the sum of those of the superbubbles. We first consider the averaged volume shocked by SN feedback per SN considering the variation of superbubbles created by a range of star clusters:

Then the effective radius is obtained as

and the shock radius by ${\cal N}_{\rm SN}$ SNe is calculated as

Here, ${\cal N}_{\rm SN}$ is meant to be the number of SNe occurring in the stellar population of a star particle. Since Eq. (34) gives the averaged superbubble volume per SN, the effective volume for $\mathcal{N}_{\rm SN}$ SNe for a star particle would be $\mathcal{N}_{\rm SN}\hat{V}$ and the effective radius would be $\mathcal{N}_{\rm SN}^{1/3}\hat{R}$ as given in Eq. (36).

Another way to determine the shock radius is to use the radius of a single superbubble described in Eq. (26) at

assuming that the percolated superbubbles from multiple star clusters can be approximated well as a single superbubble. Then, inserting Eq. (37) into Eq. (26) and using Eq. (30) with $E_{51}=1$ yields the shock radius

Similarly to the above, one can interpret Eq. (36) as the shock radius of a single superbubble formed by ${\cal N}_{\rm SN}$ SNe at the age of $t=4.8\,{\rm Myr}\,{\cal N}_{\rm SN}^{0.38}\,n_{0}^{0.22}$, which can be obtained by equating Eq. (26) and (36). In other words, the above two different methods of computing the effective shock radius can be interpreted as single superbubble radius at different times.

In the following, we use Eq. (36) together with Eq. (10) and Eq. (33) for our SN feedback model described in the next section, but the choice of either Eq. (36) or (38) is unlikely to affect our final results.

A schematic description of the ‘Spherical superbubble model’ developed in this paper is illustrated in Fig. 10. When the SN event occurs, we first calculate the density and metallicity at the SN site. An iterative solver is used to find a smoothing length for the stellar particle that satisfies

where $r_{ij}$ is the distance from the $i$-th stellar particle to the $j$-th gas particle, $N_{\rm ngb}$ is the number of neighboring SPH particles, and $W$ is the kernel function adopted in the SPH simulation. We then compute the shock radius $R_{\rm shock}$ using local density and metallicity, using equation (36). We search for the gas particles within the shock radius and project them from the position of the star onto a sphere centered at the star with radius $R_{\rm shock}$. Then, we construct a Voronoi polyhedron using STRIPACK⁶⁶6https://people.sc.fsu.edu/~jburkardt/f_src/stripack/stripack.html (Renka, 1997), which is an algorithm that constructs a Voronoi diagram of a set of points on the surface of a sphere. After constructing the Voronoi polyhedron on the spherical surface, we calculate $\Omega$, which is the solid angle of the corresponding face on the Voronoi polyhedron from the star. We distribute physical quantities from the SN to neighboring gas particles weighted by $\Omega$. The mass and metal deposited on the $i$-th gas particle are

where $m_{\mathrm{SN}}$ and $m_{Z,\mathrm{SN}}$ are the mass and metal inputs from feedback and $\Omega_{i}$ is the solid angle of the corresponding face on the Voronoi polyhedron from the star. If the number of gas particles inside the shock radius is less than four, the Voronoi polyhedron cannot be constructed. In that case, we search for at least two nearest gas particles and equally assign mass and metals to them. Note that our method is similar to, but different in detail from, that of Hopkins et al. (2018).

Using the superbubble momentum $p_{\rm SB}$ computed in Eq. (33), we deposit the following momentum on the $i$-th gas particle:

where $\bm{n}_{i}$ is the normal vector of the face on the Voronoi polyhedron. When the number of neighboring gas particles falls below four (which prevents the construction of a Voronoi polyhedron), we inject the same amount of momentum and determine $\bm{n}_{i}$ so that the total momentum is conserved. The total momentum of the surrounding gas should be conserved before and after the SN event. For this, we compute the total momentum input as:

and modify momentum input to $i$-th gas particle to

The kinetic energy input by momentum kick may exceed the SN energy input due to particle distribution. Thus, we limit the momentum input to each neighboring particle based on the Sedov–Taylor solution. The resulting momentum input is:

where $m_{i}$ is the mass of the $i$-th gas particle and $\Delta E_{\mathrm{kin},i}$ is the solid-angle-weighted kinetic energy from SN feedback. $\Delta E_{\mathrm{kin},i}$ is given as

where $\epsilon_{\rm kin}$ is a fraction of the SN energy deposited as kinetic energy, and we adopt $\epsilon_{\rm kin}=0.3$ as the default value (S19). Equation (45) essentially corresponds to equation (A1) in Kimm & Cen (2014) and equation (32) in Hopkins et al. (2018).

The momentum input above is calculated with respect to the frame moving with the star particle. To ensure exact conservation, we require a term accounting for the relative motion between the gas and the star (Hopkins et al., 2018). Finally, the momentum input boosted back to the simulation frame is

where $\bm{v}_{\rm star}$ is the star velocity.

In summary, we first compute $\Delta\bm{p}_{i}$ for the momentum that we want to assign the neighboring gas particles with. Then $\Delta\bm{p}_{i}^{\prime}$ makes it isotropic and we ensure energy conservation by $\Delta\bm{p}^{\prime\prime}_{i}$. Finally, $\Delta\bm{p}^{\prime\prime\prime}_{i}$ takes care of the motion of the originating stellar particle. By giving $\Delta\bm{p}^{\prime\prime\prime}_{i}$, the momentum feedback is basically guaranteed to be isotropic, energy-conserving, and momentum-conserving.

To be more specific, it is possible that momentum conservation could be broken when the second term in the RHS of Eq. (45) is chosen. However, such a situation does not arise very often because we do not have sufficient to resolve the Sedov–Taylor phase. In fact, we have checked our isolated galaxy simulations performed in Section 4, the second term was chosen only $3.4\times 10^{-5}$ % of the cases for the Fiducial run, and $1.1\times 10^{-4}\,\%$ for the High-reso run.

Several groups have studied the SN-driven outflow using high-resolution small-box simulations. Hu (2019) investigated SN-driven outflow of a dwarf galaxy. They showed the entropy $S\equiv{k_{\rm B}}Tn^{1-\gamma}$ of the hot outflow to be $10^{8}$ – $10^{9}\,{k_{\rm B}}\,{\rm K}\,\mathrm{cm}^{2}$, and it is almost constant after being launched from the galaxy. Although their result is on a dwarf galaxy, the energy loading factor and specific energy of the hot outflow by Hu (2019) are similar to those of a Milky-way-mass galaxy (Li et al., 2017; Armillotta et al., 2019; Kim & Ostriker, 2018; Kim et al., 2020a), according to Li & Bryan (2020).

It is suggested that the hot outflow is driven by buoyancy and its entropy is its fundamental physical quantity; if the entropy of the hot bubble is higher than that of the surrounding CGM, the hot bubble becomes buoyant and drives the outflow (Bower et al., 2017). Keller et al. (2020) analytically calculated the entropy of superbubbles and the virialized Milky-way-mass halos to be about $10^{8}\,{k_{\rm B}}\,{\rm K}\,\mathrm{cm}^{2}$. The buoyancy-driven outflow framework is supported by the MUGS2 simulations (Keller et al., 2015, 2016), and the framework may explain the ineffectiveness of SN feedback in halos more massive than $10^{12}\,{\rm M_{\odot}}$ (Keller et al., 2020).

In this work, we update the stochastic thermal feedback model (Dalla Vecchia & Schaye, 2012) by using the entropy of hot outflow $S_{\rm OF}$ as a free parameter, setting $S_{\rm OF}=10^{8}{k_{\rm B}}\,{\rm K}\,\mathrm{cm}^{2}$ as the default value. When SN feedback occurs, thermal energy is stochastically injected to heat neighboring particles to target entropy, $S_{\rm OF}$. The thermal energy required to heat the $i$-th gas particle is

where $n_{i}$ is the number density of the $i$-th gas particle. The probability of injecting thermal energy is the ratio between the solid-angle-weighted thermal energy from SN feedback (Section 3.3),

and the required thermal energy:

where $\epsilon_{\rm th}=1-\epsilon_{\rm kin}$ is a fraction of the SN energy deposited as thermal energy, and we adopt $\epsilon_{\rm th}=0.7$ as the default value (S19). We draw a random number $0<\theta_{i}<1$ for each gas particle, and inject thermal energy $\Delta E_{{\rm req},i}$ to $i$-th gas particle if $\theta_{i}<P_{i}$. When $P_{i}>1$, we simply inject a thermal energy of $\Delta E_{\mathrm{th},i}$ to the $i$-th gas particle.

The original stochastic thermal feedback model by Dalla Vecchia & Schaye (2012) uses temperature increase $\Delta T$ as a free parameter and stochastically heats gas that magnitude. Their default value $\Delta T=10^{7.5}\,{\rm K}$ was set for a numerical reason; the expectation value for the number of heated gas particles per SN is about 1. If this value is much smaller than 1, most SN feedback events do not inject energy into their surroundings, leading to poor sampling of the SN feedback cycle (Schaye et al., 2015). Although they provide a good reason to use $\Delta T$ as a parameter for their stochastic feedback model, the outflow properties depend on $\Delta T$ (Dalla Vecchia & Schaye, 2012) and a better parameter choice is desired. In this work, we use outflow entropy as a free parameter, motivated by the high-resolution simulations (Hu, 2019) and buoyancy-driven outflow framework (Bower et al., 2017; Keller et al., 2020).

We use the GADGET3-Osaka cosmological smoothed particle hydrodynamics (SPH) code (Aoyama et al., 2017; Shimizu et al., 2019), which is a modified version of GADGET-3 (originally described in Springel 2005, as GADGET-2⁷⁷7https://wwwmpa.mpa-garching.mpg.de/gadget/). We solve the SPH equation of motion, following the entropy-conserving, density-independent SPH formulation (Hopkins, 2013; Saitoh & Makino, 2013):

where $m_{i}$, $\bm{v}_{i}$, $S_{i}$, $h_{\mathrm{sml},i}$, and $\bar{P}_{i}$ denote the mass, velocity, entropy (defined as $S\equiv\bar{P}/\rho^{\gamma}$), smoothing length, and smoothed pressure defined as

of the $i$-th gas particle, respectively. In this formulation, the thermal energy of a particle is computed from its entropy and smoothed pressure. When thermal feedback injects energy, updating entropy assuming fixed pressure leads to the wrong result because the smoothed pressure field itself depends on the entropy. Therefore, we use an iterative method to calculate entropy change by energy injection from feedback and energy dissipation by radiative cooling (Schaye et al., 2015; Borrow et al., 2021). The self-gravity of SPH and collisionless particles are also considered. We adopt the quintic B-spline kernel (Schoenberg, 1946; Morris, 1996), and set the number of neighboring particles for each SPH particle to $N_{\rm ngb}=128\pm 8$. Our code includes the time-step limiter (Saitoh & Makino, 2009; Durier & Dalla Vecchia, 2012).

Radiative cooling is calculated using the grackle-3 chemistry and cooling library⁸⁸8https://grackle.readthedocs.io
The default solar metallicity in grackle-3 is $Z_{\odot}=0.0134$. This value is smaller than the value of $Z_{\odot}$ that was used in Eq. (10) by a factor of 1.45. This difference will affect the calculation of terminal momentum in Eq. (16), however the power index on $\Lambda_{6,-22}$ term is small so that the impact on the value of $p_{\rm sf,m}$ is negligible. (Smith et al., 2017), which solves non-equilibrium primordial chemistry and cooling for the H, D, and He species, including molecular H₂ and HD. The library also includes tabulated rates of metal cooling calculated with the photoionization code Cloudy (Ferland et al., 2013) and photoheating and photoionization from the ultraviolet background (UVB) radiation, and we adopt the UVB value at $z=0$ by Haardt & Madau (2012). We applied a nonthermal Jeans pressure floor that forces the local Jeans length to be resolved to avoid artificial numerical fragmentation (Hopkins et al., 2011; Kim et al., 2016):

where $\gamma=5/3$ is the adiabatic index, $N_{\rm Jeans}=4$ is the Jeans number adopted from Truelove et al. (1997), $G$ is the gravitational constant, and $\rho_{\rm gas}$ is the gas density. $\Delta x$ is chosen from the larger one of either the gravitational softening length of an SPH particle or the spatial resolution of hydrodynamics $(m_{\rm gas}/\rho_{\rm gas})^{1/3}$, where $m_{\rm gas}$ is the mass of the gas particle.

We used an initial condition taken from the AGORA project⁹⁹9http://www.AGORAsimulations.org (Kim et al., 2016). The galaxy has properties characteristic of Milky Way-mass galaxies at redshift $z\sim 1$. The galaxy is composed of the following components: a dark matter halo with $M_{\rm DM}=1.25\times 10^{12}\,{\rm M_{\odot}}$, a stellar disk with $M_{\rm disc}=4.30\times 10^{9}\,{\rm M_{\odot}}$, a stellar bulge with $M_{\rm bulge}=3.44\times 10^{10}\,{\rm M_{\odot}}$, and a gas disk with $M_{\rm gas}=8.59\times 10^{9}\,{\rm M_{\odot}}$. The total mass of the galaxy is $1.3\times 10^{12}\,{\rm M_{\odot}}$. In the fiducial run, we employed $10^{5}$ dark matter particles, $10^{5}$ gas (SPH) particles, and $10^{5}$ and $1.25\times 10^{4}$ collisionless particles representing the stars in the disk and the bulge, respectively. We also added a hot gaseous halo following Shin et al. (2021). We randomly sampled $4\times 10^{4}$ dark matter particles and added gas particles with the same mass as that of originally existing gas particles at the same position as those sampled dark matter particles. We adopt a fixed gravitational softening length of $\epsilon_{\rm grav}$ = 80 pc. We allowed the minimum gas smoothing length to reach 10 percent of the spline size of gravitational softening $2.8\epsilon_{\rm grav}$.¹⁰¹⁰10We use the definitions of the smoothing length and the gravitational softening length in GADGET-2 code (Springel, 2005) in this paper; the gravitational softening length is equivalent to the Plummer softening length, while the smoothing length is the kernel size beyond which kernel value vanishes (see also Appendix C in Kim et al., 2016). We define the SPH spatial resolution as the smoothing length widely used in literature, $\Delta x=\eta(m/\rho)^{1/3}$ (see e.g., Rosswog, 2009, for review). The parameter $\eta$ is usually set to $\eta\sim 1.3$, but we set it to $\eta=1$ for simplicity. In our GADGET3-Osaka simulation, we adopt the quintic spline kernel and $N_{\rm ngb}=128$. If we consider the quintic spline kernel to truncate at $3h$ as in Price (2012) (see also Dehnen & Aly, 2012, for another definition of the smoothing length in terms of the smoothing kernel), we set to $\eta=1.04$ in our simulation, effectively. We first evolve the system to 0.5 Gyr adiabatically for relaxation to set up the initial condition, and then evolve it to 1 Gyr with the sub-grid physics including cooling, UVB heating, star formation, and stellar feedback.