Hot phase generation by supernovae: resolution, chemistry and thermal conduction

TL;DR

This study investigates the resolution requirements for accurately simulating supernova-driven hot gas, chemistry, and momentum in the interstellar medium, emphasizing the need for resolutions finer than 10 M$_{\ m \odot}$ in galaxy formation models.

Contribution

It provides detailed resolution benchmarks for supernova simulations including chemistry and thermal conduction, highlighting the importance of high resolution for realistic ISM modeling.

Findings

Hot phase and chemistry require higher resolution than momentum.

Resolutions finer than 10 M$_{\odot}$ are necessary for accurate chemistry and hot gas modeling.

Thermal conduction reduces hot phase mass by about 0.2 dex.

Abstract

Supernovae (SN) generate hot gas in the interstellar medium (ISM), help setting the ISM structure and support the driving of outflows. It is important to resolve the hot gas generation for galaxy formation simulations at solar mass and sub-parsec resolution which realise individual supernova (SN) explosions with ambient densities varying by several orders of magnitude in a realistic multi-phase ISM. We test resolution requirements by simulating SN blast waves at three metallicities ($Z = 0.01, 0.1$ and $1 Z_{\odot}$), six densities and their respective equilibrium chemical compositions ($n=0.001$ cm$^{-3}$ - $100$ cm$^{-3}$), and four mass resolutions ($0.1$ - $100$ M$_{\odot}$), in three dimensions. We include non-equilibrium cooling and chemistry, a homogenous interstellar radiation field, and shielding with a modern pressure-energy smoothed particle hydrodynamics (SPH) method including isotropic thermal conduction and a meshless-finite-mass (MFM) solver. We find stronger resolution requirements for chemistry and hot phase generation than for momentum generation. While at $10$ M$_{\odot}$ the radial momenta at the end of the Sedov phase start converging, the hot phase generation and chemistry require higher resolutions to represent the neutral to ionised hydrogen fraction at the end of the Sedov phase correctly. Thermal conduction typically reduces the hot phase by $0.2$ dex and has little impact on the chemical composition. In general, our $1$, and $0.1$ M$_{\odot}$ results agree well with previous numerical and analytic estimates. We conclude that for the thermal energy injection SN model presented here resolutions higher than $10$ M$_{\odot}$ are required to model the chemistry, momentum and hot phase generation in a multi-phase ISM.

Tables

Table 1. Table 1: Simulation parameters for supernovae in homogeneous ambient media

Simulation Name	ambient n [cm-3]	ambient T [K]	metallicity [Z⊙]	chemistry	mgas [M⊙]	solver	SN-feedback
							resolved
Cold Molecular Medium
CMM-L1	$100$	$10$	$0.1$	$100$% H2	$100$	MFMH	\textcolorred✗
CMM-L2	$100$	$10$	$0.1$	$100$% H2	$10$	MFM	$\bullet$
CMM-L3	$100$	$10$	$0.1$	$100$% H2	$1$	MFM	\textcolorgreen✓
CMM-L4	$100$	$10$	$0.1$	$100$% H2	$0.1$	MFM, SPH	\textcolorgreen✓
CMM-cond-L4	$100$	$10$	$0.1$	$100$% H2	$0.1$	SPH	\textcolorgreen✓
CMM-Z1-L1	$100$	$10$	$1$	$100$% H2	$100$	MFM	\textcolorred✗
CMM-Z1-L2	$100$	$10$	$1$	$100$% H2	$10$	MFM	$\bullet$
CMM-Z1-L3	$100$	$10$	$1$	$100$% H2	$1$	MFM	\textcolorgreen✓
CMM-Z1-L4	$100$	$10$	$1$	$100$% H2	$0.1$	MFM	\textcolorgreen✓
CMM-Z001-L1	$100$	$10$	$0.01$	$100$% H2	$100$	MFM	\textcolorred✗
CMM-Z001-L2	$100$	$10$	$0.01$	$100$% H2	$10$	MFM	$\bullet$
CMM-Z001-L3	$100$	$10$	$0.01$	$100$% H2	$1$	MFM	\textcolorgreen✓
CMM-Z001-L4	$100$	$10$	$0.01$	$100$% H2	$0.1$	MFM	\textcolorgreen✓
Cold Neutral Medium
CNM-n10-L1	$10$	$100$	$0.1$	$100$% H	$100$	MFM	\textcolorred✗
CNM-n10-L2	$10$	$100$	$0.1$	$100$% H	$10$	MFM	$\bullet$
CNM-n10-L3	$10$	$100$	$0.1$	$100$% H	$1.0$	MFM	\textcolorgreen✓
CNM-n10-L4	$10$	$100$	$0.1$	$100$% H	$0.1$	MFM, SPH	\textcolorgreen✓
CNM-n10-cond-L4	$10$	$100$	$0.1$	$100$% H	$0.1$	SPH	\textcolorgreen✓
CNM-L1	$1$	$100$	$0.1$	$100$% H	$100$	MFM, SPH	\textcolorred✗
CNM-L2	$1$	$100$	$0.1$	$100$% H	$10$	MFM, SPH	$\bullet$
CNM-L3	$1$	$100$	$0.1$	$100$% H	$1.0$	MFM, SPH	\textcolorgreen✓
CNM-L4	$1$	$100$	$0.1$	$100$% H	$0.1$	MFM, SPH	\textcolorgreen✓
CNM-cond-L4	$1$	$100$	$0.1$	$100$% H	$0.1$	SPH	\textcolorgreen✓
CNM-Z1-L1	$1$	$100$	$1$	$100$% H	$100$	MFM, SPH	\textcolorred✗
CNM-Z1-L2	$1$	$100$	$1$	$100$% H	$10$	MFM, SPH	$\bullet$
CNM-Z1-L3	$1$	$100$	$1$	$100$% H	$1.0$	MFM, SPH	\textcolorgreen✓
CNM-Z1-L4	$1$	$100$	$1$	$100$% H	$0.1$	MFM, SPH	\textcolorgreen✓
CNM-Z001-L1	$1$	$100$	$0.01$	$100$% H	$100$	MFM, SPH	\textcolorred✗
CNM-Z001-L2	$1$	$100$	$0.01$	$100$% H	$10$	MFM, SPH	$\bullet$
CNM-Z001-L3	$1$	$100$	$0.01$	$100$% H	$1.0$	MFM, SPH	\textcolorgreen✓
CNM-Z001-L4	$1$	$100$	$0.01$	$100$% H	$0.1$	MFM, SPH	\textcolorgreen✓
Warm Neutral Medium
WNM-L1	$0.1$	$4000$	$0.1$	$100$% H	$100$	MFM	\textcolorred✗
WNM-L2	$0.1$	$4000$	$0.1$	$100$% H	$10$	MFM	\textcolorgreen✓
WNM-L3	$0.1$	$4000$	$0.1$	$100$% H	$1.0$	MFM	\textcolorgreen✓
WNM-L4	$0.1$	$4000$	$0.1$	$100$% H	$0.1$	MFM, SPH	\textcolorgreen✓
WNM-cond-L4	$0.1$	$4000$	$0.1$	$100$% H	$0.1$	SPH	\textcolorgreen✓
WNM-Z1-L1	$0.1$	$4000$	$1$	$100$% H	$100$	MFM, SPH	\textcolorred✗
WNM-Z1-L2	$0.1$	$4000$	$1$	$100$% H	$10$	MFM	\textcolorgreen✓
WNM-Z1-L3	$0.1$	$4000$	$1$	$100$% H	$1.0$	MFM	\textcolorgreen✓
WNM-Z1-L4	$0.1$	$4000$	$1$	$100$% H	$0.1$	MFM	\textcolorgreen✓
WNM-Z001-L1	$0.1$	$4000$	$0.01$	$100$% H	$100$	MFM	\textcolorred✗
WNM-Z001-L2	$0.1$	$4000$	$0.01$	$100$% H	$10$	MFM	\textcolorgreen✓
WNM-Z001-L3	$0.1$	$4000$	$0.01$	$100$% H	$1.0$	MFM	\textcolorgreen✓
WNM-Z001-L4	$0.1$	$4000$	$0.01$	$100$% H	$0.1$	MFM, SPH	\textcolorgreen✓
Warm Ionised Medium
WIM-L1	$0.01$	$8000$	$0.1$	$100$% H+	$100$	MFM	\textcolorred✗
WIM-L2	$0.01$	$8000$	$0.1$	$100$% H+	$10$	MFM	\textcolorgreen✓
WIM-L3	$0.01$	$8000$	$0.1$	$100$% H+	$1.0$	MFM	\textcolorgreen✓
WIM-L4	$0.01$	$8000$	$0.1$	$100$% H+	$0.1$	MFM, SPH	\textcolorgreen✓
WIM-cond-L4	$0.01$	$8000$	$0.1$	$100$% H+	$0.1$	SPH	\textcolorgreen✓

Simulation Name	ambient n [cm-3]	ambient T [K]	metallicity [Z⊙]	chemistry	mgas [M⊙]	solver	SN-feedback
							resolved
Warm Ionised Medium (continued)
WIM-Z1-L1	$0.01$	$8000$	$1$	$100$% H+	$100$	MFM	\textcolorred✗
WIM-Z1-L2	$0.01$	$8000$	$1$	$100$% H+	$10$	MFM	\textcolorgreen✓
WIM-Z1-L3	$0.01$	$8000$	$1$	$100$% H+	$1.0$	MFM	\textcolorgreen✓
WIM-Z1-L4	$0.01$	$8000$	$1$	$100$% H+	$0.1$	MFM	\textcolorgreen✓
WIM-Z001-L1	$0.01$	$8000$	$0.01$	$100$% H+	$100$	MFM	\textcolorred✗
WIM-Z001-L2	$0.01$	$8000$	$0.01$	$100$% H+	$10$	MFM	\textcolorgreen✓
WIM-Z001-L3	$0.01$	$8000$	$0.01$	$100$% H+	$1.0$	MFM	\textcolorgreen✓
WIM-Z001-L4	$0.01$	$8000$	$0.01$	$100$% H+	$0.1$	MFM	\textcolorgreen✓
WIM-n0001-L1	$0.001$	$10000$	$0.1$	$100$% H+	$100$	MFM, SPH	\textcolorred✗
WIM-n0001-L2	$0.001$	$10000$	$0.1$	$100$% H+	$10$	MFM	\textcolorgreen✓
WIM-n0001-L3	$0.001$	$10000$	$0.1$	$100$% H+	$1.0$	MFM	\textcolorgreen✓
WIM-n0001-L4	$0.001$	$10000$	$0.1$	$100$% H+	$0.1$	MFM	\textcolorgreen✓
WIM-n0001-cond-L4	$0.001$	$10000$	$0.1$	$100$% H+	$0.1$	SPH	\textcolorgreen✓
WIM-Z1-L4	$0.001$	$10000$	$1$	$100$% H+	$0.1$	MFM, SPH	\textcolorgreen✓
WIM-Z001-L4	$0.001$	$10000$	$0.01$	$100$% H+	$0.1$	MFM, SPH	\textcolorgreen✓

Equations

$$\mathbf{j}=-\kappa \nabla T,$$

$$\rho \frac{du}{dt}=-\nabla \text{j},$$

$$\frac{du}{dt}=\frac{1}{\rho }\nabla \cdot (\kappa \nabla T).$$

$${\kappa }_{\mathrm{sp}}=1.31n_e{\lambda }_e{k}_B{\left(\frac{k_B{T}_e}{m_e}\right)}^{1/2},$$

$$n_e{\lambda }_e=\frac{3^{3/2}(k_B{T}_e{)}^2}{4{\pi }^{1/2}{e}^4\ln \Lambda }.$$

$${\kappa }_{\mathrm{sp}}=8.2\cdot 10^{20}{\left(\frac{k_{\mathrm{B}}T}{10\text{keV}}\right)}^{5/2}\frac{\text{erg}}{\text{cm s keV}}.$$

$$j_{\mathrm{sat}}\approx 0.4n_e{k}_{B}T{\left(\frac{2k_{B}T}{\pi m_e}\right)}^{1/2}.$$

$${\kappa }_{\mathrm{eff}}=\frac{{\kappa }_{\mathrm{sp}}}{1+4.2{\lambda }_e/l_T},$$

$$m_{\mathrm{resolved}}=N_{\mathrm{ngb}}{m}_{\mathrm{gas}}{\left(\frac{h}{H}\right)}^3,$$

$$t_{\mathrm{st}}(Z=0.01Z_{\odot })=5.3\cdot 10^4\mathrm{yr}\cdot n^{-0.50},$$

$$t_{\mathrm{st}}(Z=0.1Z_{\odot })=4.2\cdot 10^4\mathrm{yr}\cdot n^{-0.50},$$

$$t_{\mathrm{st}}(Z=1Z_{\odot })=2.9\cdot 10^4\mathrm{yr}\cdot n^{-0.53}.$$

$$p_{\mathrm{st}}(Z=0.01)=18.0\cdot p_0\cdot n^{-0.16},$$

$$p_{\mathrm{st}}(Z=0.1)=16.4\cdot p_0\cdot n^{-0.18},$$

$$p_{\mathrm{st}}(Z=1)=12.0\cdot p_0\cdot n^{-0.17}.$$

$$r_{\mathrm{st}}(Z=0.01Z_{\odot })=18.0\mathrm{pc}{n}^{-0.53}$$

$$r_{\mathrm{st}}(Z=0.1Z_{\odot })=16.4\mathrm{pc}{n}^{-0.52}$$

$$r_{\mathrm{st}}(Z=1Z_{\odot })=13.15\mathrm{pc}{n}^{-0.52}$$

$${ϵ}_{\mathrm{hot}}=S\frac{4\pi }{3}{r}_{st}^3{t}_{st},$$

$${ϵ}_{\mathrm{hot}}(Z=0.01Z_{\odot })=S\cdot 1.310^{-6}{\mathrm{kpc}}^3{\mathrm{Myr}}^{-1}{n}^{-2.1},$$

$${ϵ}_{\mathrm{hot}}(Z=0.1Z_{\odot })=S\cdot 9.7410^{-7}{\mathrm{kpc}}^3{\mathrm{Myr}}^{-1}{n}^{-2.1},$$

$${ϵ}_{\mathrm{hot}}(Z=1Z_{\odot })=S\cdot 6.6210^{-7}{\mathrm{kpc}}^3{\mathrm{Myr}}^{-1}{n}^{-2.1}.$$

$$M_{\mathrm{hot}}(Z=0.01Z_{\odot })=2570\cdot M_{\odot }\cdot n^{-0.29},$$

$$M_{\mathrm{hot}}(Z=0.1Z_{\odot })=2398\cdot M_{\odot }\cdot n^{-0.29},$$

$$M_{\mathrm{hot}}(Z=1Z_{\odot })=1584\cdot M_{\odot }\cdot n^{-0.28},$$

$$t_{\mathrm{st},\mathrm{sph}}=3.8\cdot 10^4\mathrm{yr}\cdot n^{-0.50},$$

$$t_{\mathrm{st},\mathrm{cond}}=3.6\cdot 10^4\mathrm{yr}\cdot n^{-0.46}.$$

$$p_{\mathrm{st},\mathrm{sph}}=14.6\cdot p_0\cdot n^{-0.18},$$

$$p_{\mathrm{st},\mathrm{cond}}=12.9\cdot p_0\cdot n^{-0.18}.$$

$$M_{\mathrm{hot},\mathrm{cond}}=1580\cdot M_{\odot }\cdot n^{-0.29}.$$

Full text

††pubyear: 2019††pagerange: Hot phase generation by supernovae in ISM simulations: resolution, chemistry and thermal conduction–A

1 Introduction

The galactic interstellar medium is shaped by the feedback of massive stars. A physical process which is of paramount importance in this picture are SN-explosions. Core-collapse-supernovae (SNe) eject gas at supersonic velocities of several thousand kilometres per second (e.g. Blondin et al., 1998; Janka et al., 2012) and drive blast waves into the ISM. The evolution of SN-remnants has been extensively investigated analytically, and can be divided in four phases (e.g. Woltjer, 1972): the phase of ’free-expansion’ which is terminated when the swept up mass equals the ejecta mass (roughly a few 100 years), the energy conserving Sedov-Taylor (ST) phase (e.g. Sedov, 1946, 1959; Taylor, 1950), the pressure driven snow plough (PDS) phase and the momentum conserving snow plough phase (MCS; Falle 1975). In the PDS-phase the pressure of the bubble interior further generates momentum until pressure equilibrium with the shell is reached and the bubble evolves the remnant in the low pressure regime (Cox, 1972; Gaffet, 1983; Cioffi et al., 1988; Cohen et al., 1998; Haid et al., 2016). In the MCS-phase the momentum of the shell is conserved and the further evolution of the shell is driven by the generated inertia of the swept up mass (Cioffi et al., 1988; Haid et al., 2016). The energy conserving ST-phase is of particular importance as it is the main momentum generating phase with an expanding hot bubble. This momentum build up due the ST-phase drives turbulence within the ISM and regulates star formation within galaxies. It can last from a few $1000$ years in dense environments to several million years in low density environments. The duration of the ST-phase and the momentum and hot phase generation during this phase is highly dependent on the cooling properties of the ambient ISM. In this phase the SN-remnant behaves adiabatically because there are (nearly) no cooling losses that can remove energy from the system.

Blast waves replenish the hot phase of the ISM (Cox & Smith, 1974; McKee & Ostriker, 1977) through shocks (e.g. Gotthelf et al., 2001), generate momentum and build up turbulence (Elmegreen & Scalo, 2004; Scalo & Elmegreen, 2004) in the warm neutral medium (WNM) of the ISM (Mac Low & Klessen, 2004). The hot phase generated by SNe fills a large fraction of the turbulent ISM volume (e.g. McKee & Ostriker, 1977; Ferrière, 1998; Könyves et al., 2007) in which smaller and cooler clouds are embedded, as observed in the local ISM (e.g. Frisch et al., 2011). SNe also enrich the ISM with metals and dust (e.g. Dwek, 1998; Indebetouw et al., 2014; Kobayashi et al., 2011; Matsuura et al., 2011), different atomic species, and form molecules in their remnants (Spyromilio et al., 1988; Grefenstette et al., 2014; Fransson et al., 2016; Kamenetzky et al., 2013). When multiple SNe occur in a low density environment, subsequent explosion add mass to the shell and heat the interior of the bubble which can lead to the formation of ’super-bubbles’ (Castor et al., 1975; Weaver et al., 1977; McCray & Kafatos, 1987; Tomisaka & Ikeuchi, 1986; Mac Low & McCray, 1988; Koo & McKee, 1992). Theses can drive strong galactic outflows which redistribute the gas on galactic scales and moves a fraction of it towards the circum galactic medium (CGM) of the galaxy. In more general terms, the feedback from SNe can regulate the cosmic baryon cycle and the star formation rate (SFR) of galaxies across cosmic time (c.f. Somerville & Davé, 2015; Naab & Ostriker, 2017).

In cosmological simulations of galaxy formation, stellar feedback leads to a better agreement with the properties of the observed stellar and gaseous components, such as the low conversion efficiency of gas into stars which are determined empirically (e.g. Moster et al., 2010; Behroozi et al., 2010; Moster et al., 2013) and disk-like galaxy morphologies (Naab & Ostriker, 2017; Somerville & Davé, 2015). However, most simulations cannot resolve the momentum and hot medium generating phases of blast waves and rely on "sub-resolution" implementations to mimic the SN impact. Simulations adopting these methods have been able to produce galaxy morphology and population properties in agreement with observations (e.g. Guedes et al., 2011; Aumer & White, 2013; Marinacci et al., 2014; Agertz et al., 2013; Wang et al., 2015; Hopkins et al., 2018), and empirical predictions (e.g. Moster et al., 2018; Behroozi et al., 2018). However, it remains unclear how well the underlying physical processes are really captured, and to what degree these results are achieved by fine-tuning model parameters.

Recently, it has become possible to simulate low mass galactic systems like dwarf galaxies at such high resolution that the impact of supernovae from individual massive stars can be captured in the relevant temperature and density regimes of the multi-phase ISM (Forbes et al., 2016; Hu et al., 2016, 2017; Emerick et al., 2019). In particular studies by (Lahén et al., 2019b; Lahén et al., 2019a) have extended previous galactic ISM studies to environments with much more extreme densities and star formation rates. This approach, however, requires a detailed understanding of the resolution requirements for accurately tracking the impact of individual blast waves in the turbulent cold, warm and hot ISM. Many high-resolution simulations of SN blast waves in ambient homogeneous media or turbulent ISM patches have been carried out to understand their evolution in detail (e.g. Ostriker & McKee, 1988; Blondin et al., 1998; Thornton et al., 1998; Draine, 2011; Kim & Ostriker, 2015; Walch et al., 2015; Gatto et al., 2015; Martizzi et al., 2015; Badjin et al., 2016; Haid et al., 2016; Ohlin et al., 2019). While some of these studies aim for detailed blast-wave evolution (e.g. Thornton et al., 1998; Badjin et al., 2016) in one or two dimensional simulations, others aim for the effects of supernova feedback in a three-dimensional approach for investigating the effects on the three-dimensional galactic ISM (e.g. Kim & Ostriker, 2015; Ohlin et al., 2019) or investigate the effect of multiple feedback events on the galactic ISM (e.g. Gatto et al., 2015).

Apart from emerging resolution requirements, it has become clear that the distribution of ambient densities of SN blast waves play an important role for how efficiently they can drive turbulence and outflows (Girichidis et al., 2016; Hu et al., 2016, 2017; Gatto et al., 2017; Naab & Ostriker, 2017; Seifried et al., 2018; Lahén et al., 2019a). While the turbulent component of the ISM can be generated by adopting some feedback prescription that injects the momentum of SNe, the galactic winds are driven by the pressure that is generated by the SNe within the turbulent ISM. In simulations, many factors like spatial resolution, stellar winds, radiation, SN clustering and binary or runaway stars can affect the ambient SN density distributions (Kim et al., 2017; Gentry et al., 2017; Naab & Ostriker, 2017; Peters et al., 2017; Fielding et al., 2018). However, not only numerical implementation, resolution constraints and ambient density distribution of blast waves determine their impact but also the complexity of the physical modelling. For example magnetic fields, cosmic rays (e.g. Diesing & Caprioli, 2018; Gupta et al., 2018) and thermal conduction (e.g. Keller et al., 2014; El-Badry et al., 2019) change their evolution.

In this paper we use two particle based hydrodynamical methods: modern smoothed particle hydrodynamics (SPH) and the meshless finite mass (MFM) method (Gaburov & Nitadori, 2011; Hopkins, 2015). We couple both solvers to a non-equilibrium chemical network to estimate their ability to converge on blast wave evolution in cold, warm and hot ambient media at different numerical resolutions. Furthermore we use the SPH implementation to probe the impact of thermal conduction on individual blast waves. While these studies have been carried out extensively using different grid codes in one or two dimensions, results in three dimensions are rare and detailed studies for particle codes are missing. However, as the ISM is a highly structured, turbulent and multi phase fluid it is of importance to investigate the effects of supernovae on the galactic ISM at the highest resolutions we can achieve today in galaxy formation and evolution simulations in three dimensions to work out the limitations on the physics of the ISM in galactic scale simulations. We test our results against detailed blast wave studies carried out in one or two dimensions (e.g. Blondin et al., 1998; Thornton et al., 1998; Badjin et al., 2016) and check how well we can resolve the global evolution of the Sedov-Taylor phase with momentum injection and hot phase generation under controlled boundary conditions. We also address resolutions requirement for current and next generation galaxy formation simulations without sub-grid modelling of individual supernovae.

The paper is structured as follows. In section 2 we briefly present the numerical and physical methods that are important for this study and discuss the simulation setup and the initial parameters that are important for our cooling and chemistry network. In section 3 we discuss our results for isolated SNe blast waves in a homogeneous medium (e.g. structure of the blast wave its environmental dependence). In section 4 we use the physical properties from our simulations to constrain the expectation value for a hot phase to form as a function of their environmental densities. In section 5 we discuss the effect of isotropic heat conduction on the SN-remnant evolution in our highest resolution runs. In section 6 we derive the consequences of the presented supernova-feedback scheme on the ISM in galaxy scale simulations and show results for supernova-driven ISM-patches, remnants in structured media and finally full galactic disc simulations at solar mass and sub-parsec resolution. Finally, we summarise our results in section 7.

2 Simulation method

2.1 Hydrodynamics

We run our set of simulations with our version of Gadget-3 (Springel, 2005), which includes implementations for several hydrodynamics solvers, such as pressure-energy SPH (Hu et al., 2014) and the meshless finite mass (MFM) (e.g. Gaburov & Nitadori, 2011; Hopkins, 2015). In Hu et al. (2014) it is shown that the improved pressure-energy SPH accurately captures shocks and instabilities in several idealised test problems. While MFM shares the kernel weighted density computation with SPH, the hydrodynamical flux vectors are integrated over the one-dimensional Riemann-problem defined on the surface between two reconstruction points of the fluid equations (particles). In combination with a second order reconstruction of the flux gradients (Gaburov & Nitadori, 2011; Mocz et al., 2014) it is possible to obtain a second order hydrodynamical scheme with an appropriate mathematical consistency and convergence proof as shown in Lanson & Vila (2008). We follow the slope limiting and reconstruction procedures by Hopkins (2015). We calculate the quadrature point in first order via $\mathbf{x}_{ij}=(\mathbf{x}_{i}+\mathbf{x}_{j})/2$ and solve the one dimensional Riemann-problem with a Harten-Lax-van-Leer (HLL) Riemann-solver with an approximate reconstruction of the contact wave (HLLC) following Toro et al. (1994). The time integration scheme that we adopt follows the description of Springel (2010) with a common CFL-criterion for the timestep. In strong shocks driven by supernova blast waves particles can move very fast and interact with other particles, which might be on much longer time steps and cannot react accurately. To avoid this we activate all nearby particles (within a kernel radius) and put them on the same short time step. The procedure is similar to the methods described in Saitoh & Makino (2009) and Durier & Dalla Vecchia (2012).

2.2 Thermal conduction

Some of our simulations employ a prescription for thermal conduction. We follow the implementation of Jubelgas et al. (2004) with updates for this scheme presented by Petkova & Springel (2009). We implement isotropic conduction given as a local transport process for the internal energy of a fluid tracer with the heat flux $\mathbf{j}$.

[TABLE]

where $\kappa$ is the conduction coefficient and $T$ is the temperature. The total change of energy is

[TABLE]

with the thermal energy $u$ and the density of the fluid $\rho$. The conduction equation can be obtained by inserting equation 1 into 2:

[TABLE]

However, we note that in astrophysical plasmas the conductivity is not constant but shows a dependence on the temperature. In the Spitzer limit this is given by

[TABLE]

with the electron density ne, its temperature Te and its mass me. However, $n_{e}\lambda_{e}$ is only a function of the temperature Te and the Coulomb-logarithm following

[TABLE]

In the case of a constant value for the Coulomb-logarithm the Spitzer conductivity remains strongly temperature dependent

[TABLE]

However, the behaviour changes for very low plasma densities where the scale length of the temperature gradient is of order the electron mean free path or even smaller. In this regime the heat flux saturates and no longer increases even if the temperature gradient rises. The saturation flux is given by

[TABLE]

In this case $\kappa$ is modified to an effective $\kappa_{\mathrm{eff}}$ to smoothly change between the Spitzer and the saturated regime in the following manner

[TABLE]

where lT is the length scale of the temperature gradient given as $T/|\nabla T|$. While the saturation limit is mostly relevant on galaxy-cluster scales and not in the ISM, we carried out some test simulations without the heat flux limit and even without a temperature dependence of $\kappa$ itself. For isolated events we do not find large differences and it seems to be more important to have a heat diffusion term in the first place. The reason for that is, that we cannot resolve the process of thermal conduction within the bubble or the shell respectively. However, what we can resolve is the interface between bubble and shell. Therefore, even a constant value for $\kappa$ leads to similar results (as long as it is physically motivated), as $\kappa$ at the interface of bubble and shell stays (roughly) constant.

2.3 Chemistry model

We adopt a non-equilibrium chemical model following the implementation of the ’SILCC’ project (Walch et al., 2015; Girichidis et al., 2016; Gatto et al., 2017; Peters et al., 2017). This chemical model is based on Nelson & Langer (1997), Glover & Mac Low (2007a), Glover & Mac Low (2007b) and Glover & Clark (2012). We track six individual species, H2, H*+, CO, H, C+* and O as well as free electrons and assume that silicon is present in the simulation as Si*+. Carbon is present either in single ionised from as C+* or CO. Oxygen is present as O. H2 and H*+* undergo several reactions, which are summarised in table 1 of Micic et al. (2012). The formation of H2 is a dust catalysed reaction in which a dust grain captures two hydrogen atoms and the reaction becomes possible on dust surfaces. Most important for the destruction of H2 is photo-dissociation due to the inter-stellar radiation field (ISRF), ionisation due to cosmic rays and via collisional dissociation (collisions with other species). The main channels for the formation of H*+* are collisional excitation with cosmic rays and free electrons. The destruction of H*+* is driven by recombination processes with free electrons. In the case of thermal conduction we directly assume the amount of free electrons that is predicted from solving the non-equilibrium rate equations to perform our computation of heat transfer.

2.4 Cooling and heating processes

The non-equilibrium cooling rates take into account the local density, temperature and the local chemical abundances. There are six major non-equilibrium cooling processes that are relevant for our studies, fine structure line cooling of atoms and ions (C*+, Si+, O), vibration and rotation line cooling of molecules (H2*, CO), Lyman-alpha cooling of hydrogen, collisional dissociation of H2, collisional ionisation of hydrogen and recombination of H*+* (both in gas phase and condensed on dust grains). The main sources for heating are given by photo-electric heating from dust grains and polycyclic aromatic hydrocarbonates (PAHs), cosmic ray ionisation, photo-dissociation of H2, UV-pumping of H2 and formation of H2. For high temperatures ($T>3\cdot 10^{4}$ K) we adopt the cooling function by Wiersma et al. (2009). It assumes that the ISM is in equilibrium where the cooling due to collisional ionisation is balanced by the heating of the cosmic UV-radiation (Haardt & Madau, 2001). Further, the absolute cooling rate is determined by following the model of Aumer & White (2013) and depends on the local gas density, gas temperature and includes metal line cooling from 11 species (H, He, C, N, O, Ne, Mg, Si, S, Ca and Fe).

2.5 Shielding mechanisms from the interstellar radiation field

For the hydrogen chemistry the ISRF photo dissociation rate R${}_{\mathrm{pd,H_{2}}}$ acting on H2 is affected by the shielding of dust and the self shielding of H2. Dust has the ability to absorb photons of the ISRF field and re-emit in the infrared while absorbing the energy difference. Self-shielding by molecular hydrogen is more complicated. The process normally refers to a mechanism where the photo excitation transitions become optically thick and shield each other if UV radiation is coming from a single direction. However, in the optical thick limit (high column densities) this process is prevented by the line broadening of the Lyman and Werner bands (e.g. Gnedin & Draine, 2014). In this regime the self shielding by dust becomes the more important process. We include both processes, the self-shielding by dust and molecular hydrogen. This is implemented with the TreeCol-algorithm (Clark et al., 2012; Hu et al., 2016).

2.6 Numerical resolution and performance

Although the single particles have very small masses ranging from 0.1 M*⊙* to 100 M*⊙* this is not the effective resolution of the simulations. For methods that follow an SPH-like volume distribution the relevant spatial resolution is the size of the kernel or to be more specific the radius of compact support. For the SPH runs we use a Wendland C4 kernel with $100$ neighbours while we use a cubic spline with $32$ neighbours for the MFM runs to capture the fluid properties. The effective spatial and mass resolution is higher for the runs with MFM as the radius of compact support is smaller for the cubic spline and it can operate with a smaller number of neighbours. The resulting mass resolution can be calculated via

[TABLE]

where $h$ is the smoothing length and $H$ is the radius of compact support of the respective kernel. For the cubic spline we find the relation $h=0.55H$ and for the Wendland C4 kernel $h=0.45H$. Therefore, the resolution for MFM is better by about a factor of $2$. While the MFM method requires more memory than SPH methods (a factor of $1.5$) it is overall faster in solving the equations of hydrodynamics compared to our state-of-the-art SPH solver (for the same number of neighbours).

2.7 Initialisation and naming convention

In our default runs we adopt a metallicity of $Z=0.1$ Z*⊙* with a constant the dust to gas ratio mass ratio of 0.01. The constant UV-background (Haardt & Madau, 2012) has a normalisation factor of $G_{0}=1.7$ for the interstellar radiation field, which corresponds to the value of the solar neighbourhood. For resolution we refer to the gas particle mass. The four resolution levels are $100$ M*⊙* (L1), $10$ M*⊙* (L2), $1$ M*⊙* (L3) and $0.1$ M*⊙* (L4). The SN blast waves are simulated in the cold molecular medium (CMM) with $n=100$ cm*-3*, the cold neutral medium (CNM) with $n=1$ cm*-3*, the warm neutral medium (WMN) with $n=0.1$ cm*-3* and the warm ionised medium. Some runs have lower $Z=0.01$ $Z_{\odot}$ (Z-001) and and higher $Z=1Z_{\odot}$ (Z-1) metallicities. A different number density compared to our defaults is indicated with n and the number density in units of cm*-3*. Runs with isotropic thermal conduction are indicated by cond.

3 Individual supernova blast waves

We present a set of simulations of individual SN blast waves in six different environments at equilibrium conditions of the ISM for our cooling model. In the upper panel of Figure 1 we show the temperature as a function of hydrogen number density for different values for a (constant) interstellar radiation field and the dust to gas ratio (DTGR). The coloured circles indicate the ISM conditions for WIM, WNM, CNM and CMM. The CMM represents regimes in dense molecular clouds ($n_{H}=100$cm*-3*, $T=10$ K, fully molecular). The CNM has higher temperatures and densities such that hydrogen is present in neutral atomic form ($n_{H}=1$cm*-3*, $T=100$ K, fully neutral). The WNM is warmer but still neutral ($n_{H}=0.1$cm*-3*, $T=8000$ K, fully neutral) and the WIM is fully ionised with long cooling times (n${}_{H}=0.001$cm*-3*, T$=10000$ K, fully ionised). The pressure vs. density plot in the bottom panel highlights the pressure-unstable regime between $10^{-2}$ and $10^{-1}$ cm*-3*. We initialise each supernova event by a kernel-weighted injection of $10^{51}$ erg of thermal energy into all kernel particles. For the MFM-solver we hereby use 32 particles and cubic-spline (M4) kernel. For the respective SPH runs we use a Wendland C4 kernel with 100 neighbours. Experience has shown that MFM with a cubic spline and $32$ neighbours shows the same quality of capturing fluid flows as the SPH reference runs with 100 neighbours and a Wendland C4 kernel. In the case of SPH the higher neighbour number is needed to suppress the intrinsic noise of particles methods by smoothing over the fluid field. However, the well known and investigated pairing-instability makes it than necessary to use a higher order kernel with positive Fourier transformation. The initial phase of free expansion is very short (only a few $100$ yr) and is not modelled here. The initial conditions are summarised in Table 1.

3.1 Blast wave evolution

In Figure 2 we show the radial profiles of the hydrogen number density (top left), the temperature (top right), the relative velocity of the shock (bottom left) and the thermal pressure (bottom right) for a blast wave in the CNM at the highest resolution of $0.1$ M*⊙, which implies a spatial resolution of around $1$ pc. The different colours show the system at time $t=0.04$ Myr (black) during the ST-phase, at $t=0.16$ Myr (red) after the formation of a dense cooling shell, and at $t=1.0$ Myr (blue) when the remnant starts to merge with the ISM. In the ST-phase the maximum value of the density is around $3.53$ cm-3*. In the case of the pure adiabatic solution we would expect a maximum density of $4$ cm*-3* under constraint of the Rankine-Hugoniot jump-conditions with an adiabatic index of $5/3$ which is adopted in all our simulations.After the end of the ST-phase shell remnant has formed a thin, cooling dominated shell with a maximum density of roughly $20$ times the initial ambient density. At the end of the simulation the amplitude of the density decreases again and the remnant starts to merge with the ambient ISM. In all three stages the temperature declines as a function of the radius. After shell we can observe a bump in the temperature distribution at around $10^{4}$ K which can be identified as the dense cooling shell of the remnant. This shell grows over time as more and more gas from the initially hot bubble starts to cool. The velocities within the ST-phase are high which can lead to a velocity dispersion of roughly $10$ km s*-1* that is observed within the ISM. At later stages the remnant slows down until it finally reaches the point where it merges with the ISM. Initially, the pressure is large (especially in the bubble) but decreases about an order of magnitude within the later evolutionary stages. While the results do not compare well with studies that are carried out in fewer dimensions (Thornton et al., 1998; Blondin et al., 1998; Badjin et al., 2016) in terms of peak density in the Sedov-Taylor phase and in terms of resolving fluid instabilities that are certainly important in supernova-remnant modelling, we note that we compare well to other grid code results in three dimensions (e.g. Kim & Ostriker, 2015) who find a maximum value between 3.3 and 3.5 in their respective Sedov-Taylor phase. However, we note that we focus on low metalicity environments while (Kim & Ostriker, 2015) is focusing on solar metallicity environments. This leads to differences in blast wave evolution which we discuss in section 3.3. Moreover, we note that we observe a clear drop in the pressure between bubble and shell which indicates some contribution of the thermal instability in our highest resolution runs, although we see only a drop of $0.5$ dex while higher resolution state-of-the-art simulations of detailed blast wave evolution report $1$ dex (e.g. Badjin et al., 2016) and we do not properly resolve it at the presented resolutions throughout the paper. However, as we carry out these simulations to test which supernova-remnant properties can be modelled in current high resolution galaxy formation simulations, resolving the thermal instability in the remnant is not the main goal of the paper. Nevertheless, it would be very interesting to study the detailed effect of the thermal instability on momentum and hot phase generation in higher resolution simulations with our code in future studies. In this context, we have to consider the inability of capturing the analytic predicted value of the density and the inability of fully capturing the thermal instability as a limitation that has to be noted. However, despite this limitations which can be resolved at higher resolution or in lower dimensionality, the ability to model radial momentum input and hot phase at our target resolutions is worth noting as it enables the galaxy scale studies of outflows that are self-consistently driven by the the supernova-feedback without the limitations that common subgrid-models are tied to.

To understand the details of the blast waves evolution we investigate the remnant at different mass resolutions. In Figure 3 we show the remnant in the CNM during the ST-phase for all four mass resolutions that we investigated. We note that the shock is well resolved at the highest mass resolutions for the simulations CNM-L4 and CNM-L3. While the shock is still sufficiently resolved in the adiabatic regime in CNM-L2 it is completely unresolved for CNM-L1. This is mainly due to the fact that we have a lower amount of resolution elements in the same spatial region to actually resolve the shock properly in the adiabatic regime. Therefore, to capture the shock within the ST-phase a mass resolution of smaller than $10$ solar masses is needed for the MFM-solver. We mark runs that are resolved (green check), weakly resolved (black dot) and unresolved (red crosses) in Table 1. Further, we note differences in the cooling dominated regime. At the highest mass resolutions a thin, cooling dominated shell is forming following the end of the ST-phase. The density within this supernova shell is up to $10$ times higher than what we would expect to see in a pure adiabatic simulation of the shock front. However, with lower resolution this dense cooling shell becomes less and less resolved and remains unresolved at a mass resolution of $100$ M*⊙*. We note that this becomes even worse for a different ambient configuration. If we consider the runs in the CMM the shock is less well resolved even at higher mass resolution. In the last evolution stage of the remnant the density decreases again and starts to merge with the ambient medium for all mass resolutions. We note that in these simulations we cannot capture the phase in which the SN-remnant merges with the ambient medium as the configuration of our ambient medium is chosen to be of uniform density and zero velocity. Consequently, our ambient medium is not turbulent which would be necessary to properly capture the merging phase of the shock with the ambient medium.

3.2 Resolution effects

For an accurate simulation of a supernova blast wave for applications in galaxy scale simulations it is important to capture the formation of the swept-up shell, which also determines the cooling properties of a supernova remnant (e.g. Kim & Ostriker, 2015). For example Hu et al. (2016) have shown for SPH that the shell momentum is relatively insensitive to numerical resolution. However, the shell mass and velocity are strongly affected by higher velocities and lower shell masses at high resolution. At insufficient resolution the shell becomes too massive, slower and cools too early, resulting in a too rapid termination of the ST-phase. Similar results have been found by Kim & Ostriker (2015) for Eulerian grid simulations. However, they focused on solar metallicity environments while we focus our study on sub-solar metallicities. The material in the shell has a temperature below $2\cdot 10^{4}$ K. All material above this temperature cut belongs to the hot phase of the remnant. Additionally, we set a velocity cut of 0.1 km s*-1* to select the ambient medium which is at rest form the bubble and the shell. The shell radius is defined by the position of the $10$ particles with the highest densities. The radial shell momentum is computed by summing up the individual particle momenta, and the shell velocity is the shell momentum divided by the shell mass.

In Figure 4 we show shell radius (top left), shell momentum (top middle), shell energy (top right), split in kinetic (blue) and thermal energy (red), shell velocity (bottom left), shell mass and bubble mass (bottom middle), as well as the bubble temperature (bottom right). We present the four resolution levels with $0.1$ M*⊙* (solid line), $1.0$ M*⊙* (dashed-line), $10$ M*⊙* (dash-dot-line), $100$ M*⊙* (dash-dot-dot-line) for a medium with $n=1$ cm*-3* consisting of neutral hydrogen (CNM-L4 to CNM-L1). The small black dashes on the axis indicate the end of the ST-phase and the numerical values of the respective quantity for this point in time. We note that instead of showing the pressure in shell and bubble, we show the temperature structure within the hot bubble. The pressure within the bubble can be computed at any time by multiplying the densities within shell and bubble with the related temperatures. Moreover, we note that the lines for different resolutions start at different points in time because in the lower resolution runs the shock can not be separated from the ambient density field by the time at which the respective lines start. This is due to the fact that the shock remains unresolved until the respective points in time. In the top left panel of Figure 4 before shell formation t < tsf at the end of the Sedov phase, the radius is increasing as $r\propto t^{2/5}$ (Sedov, 1946; Taylor, 1950) (green line). At shell formation (t${}_{\mathrm{sf}}=45000$ yr) the slope changes $r\propto t^{2/7}$ for the momentum conserving phase. We define our Sedov-Taylor time tst as the (final) time when the thermal energy drops below the analytical value of $72$ per cent. The three highest resolution levels capture the shell radius at formation and its evolution within $10$ per cent. The shell momentum (top middle panel of Figure 4) is relatively insensitive to resolution and well captured at the three highest resolution levels. Further, our results are in good agreement with results from Kim & Ostriker (2015) and Hu et al. (2016). The shell velocity and mass, however, depend strongly on resolution (bottom left and middle panel of Figure 4). With higher resolution the remnant velocity increases and the remnant mass decreases as the shock is better resolved. However, we note that a part of this trend is related to our feedback scheme. We further distinguish between the total mass of the remnant (black), the mass in the cold shell (blue) and the mass in the hot bubble (red). With decreasing resolution the mass within the bubble is lower apart from the lowest resolution run where the initially heated particles generate a too high hot mass. However, as the structure of the shell is so poorly resolved in these runs those particles cool very inefficiently and the hot mass is over predicted. For lower mass resolution we inject initially in more mass, leading to more initial swept-up material. This has consequences for the shell momentum at lower resolution. In the $10$ M*⊙* run the shell momentum is too high because of this reason, which renders the shell momentum a bad tracer to determine whether the ST-phase is resolved. The evolution of shell and bubble temperature is shown in the lower right panel. At the two highest resolution levels the bubble temperatures are resolved to within $10$ per cent. At the lowest mass resolution the temperature structure of the hot bubble is unresolved.

3.3 Blastwave evolution past the ST-phase

While we focus on the evolution of the blastwaves within the ST-phase to investigate how momentum and hot mass are generated in this evolutionary state we briefly want to discuss the transition from the ST-phase to the PDS and finally the MCS phase. However, we note two major differences to other studies (e.g. Kim & Ostriker, 2015). Our main focus is on sub-solar metallicity environments, as they can be found in present day dwarf-galaxies or at higher redshifts. Moreover, our cooling-function is non static and we excluded the effect of cosmic-ray ionisation to study SN-feedback in high density and very low temperature environments. Thus, we find three times higher terminal momenta compared to solar metallicity studies with higher equilibrium temperatures in the ambient ISM. First, metal-line cooling is less efficient in lower metallicity environments. This leads to an increased momentum boost past the ST-phase as the bubble cools less efficiently. Second, as we exclude cosmic ray ionisation th equilibrium temperature is lower and we can gain a little momentum boost after the PDS-phase as the shell is pushing into the ISM until the ambient ISM reaches pressure-equilibrium with the shell. Thus our momentum input gives an upper limit for the momentum boost which will generally be lower in practice as we will discuss in 6. In Figure 5 we show two examples for the CMM (blue) and the WNM (black) without (dashed) and with (solid) cosmic ray ionisation heating of the ISM at a rate of $3\cdot 10^{-17}$ s*-1*. We see a decrease of the final momentum by roughly 20 per cent if we switch on cosmic ray ionisation heating as the ambient pressure is higher and the shell cannot generate momentum anymore while pushing into the ambient ISM. Further, the exclusion of the cosmic-ray ionisation rate has an effect on the chemistry. While the molecular hydrogen is not effected (less than the per cent level), the cosmic-ray ionisation leads to higher fractions of ionised hydrogen by a roughly $10$ per cent and subsequently reduces the fraction neutral hydrogen. There is no effect on the formation of CO. Finally, we note that at the two highest mass resolutions we start to observe the fragmentation of the shell by Ryu-Vishniac type of instabilities. While we do not claim that these instabilities are fully resolved at the presented resolution, we note that they can further affect the cooling-rates that we calculate on the fly in our simulations. A detailed modelling of the cooling rate in shell and bubble plays a crucial role in modelling a supernova-event in a galactic context. The Ryu-Vishniac instabilities (even though our results are not fully converged yet) only appear in more than one dimensional simulations and thus it is important to model the SN-feedback event in three spatial dimensions as they can influence the detailed cooling rates in our chemical model. We show this for our highest resolution simulation in a time evolution of four different snapshots, in the CMM-medium in density of 100 cm3 in Figure 6. We can see clear fragmentation of the remnant at very late times which is potentially caused by the Ryu-Vishniac instabilities. However, at this resolution the results could be also influenced by amplification of initial numerical noise of our particle distribution. We carefully checked the results by applying higher time resolution (i.e. decreasing the Courant number by a factor of 4) and find similar results. Moreover, we carried out a higher resolution simulation with a particle mass resolution of $10^{-4}$ M*⊙* and find similar structures.

3.4 Environmental dependence

The evolution of SN remnants is strongly dependent on the environment. The cooling-time scales as $n^{2}$ which has consequences for remnants in a high density environment as well as remnants in a low density environment and leads to various different shell formation times. In Figure 8 we show the ambient densities for all our SN-remnants at the highest resolution (L4) as a function of the shell formation time tsf. With increasing ambient densities the shell formation times become shorter and can last from a few $1000$ years at $100$ cm*-3* to roughly $1.5$ Myr at $0.001$ cm*-3*. As we simulate in three different metallicity regimes we show our best fits for the shell formation time individually for each metallicity as the dashed red line ($Z=1Z_{\odot}$), the solid red line ($Z=0.1Z_{\odot}$) and the dashed-dotted red line ($Z=0.01Z_{\odot}$). For $Z=0.1Z_{\odot}$ we additionally show the data from our six highest resolution runs to visualise the data with which we obtained the red solid line. The best fit relations as functions of the metallicity are given by:

[TABLE]

[TABLE]

[TABLE]

We compare these shell formation times to previous numerical and analytical results by Blondin et al. (1998), Cioffi et al. (1988), Franco et al. (1994), Petruk (2006) and Haid et al. (2016) and find good agreement with the recent results of Haid et al. (2016) and Petruk (2006).

In Figure 7 we show the blast wave evolution for all environmental densities at the highest resolution level (L4) corresponding to $0.1$ M*⊙* for different properties of the individual SN-remnants for an ambient metallicity of $Z=0.1Z_{\odot}$. With increasing ambient number density the size of the remnants at the end of the ST-phase becomes shorter and can vary over more than an order of magnitude. The momentum at the end of the ST-phase is the largest for remnants in low density environments. For all environments apart from the highest ambient density we find a trend of an increasing kinetic energy in the PDS-phase of the remnant while the thermal energy is radiated away quickly after the end of the ST-phase. The shell velocities are the highest in the low density media where the shell and bubble masses are the highest. Further, we find that the bubble temperature is higher for higher mass densities after $10$ shell formation times, while the shell temperature significantly drops due to the shorter cooling times. Finally, we note that we can resolve the ST-phase for all ambient densities at our highest mass resolution. At lower mass resolution convergence can still be obtained for the three highest mass resolutions up to an ambient density of 10 cm*-3*. For higher ambient media a mass resolution of $1$ M*⊙* is required to resolve the ST-phase. We show the results for a lower density environment and our highest ambient density environment in Appendix A.

The shell formation time terminates the radial momentum generation phase (see Haid et al. 2016 for a detailed discussion of the momentum generating phases) and the remnant transits to the momentum conserving snow plough phase. However, we note that given our cooling function and initial conditions our simulations indicate that the momentum is further increasing past pressure balance of bubble an shell is reached. However, as we have relatively low ambient temperatures due to our cooling function that does not contain the contribution of cosmic ray heating. While this provides relative extreme conditions with high cooling rates to determine how well the code can handle blast waves in rapidly cooling environments we lack comparability towards other studies like Kim & Ostriker (2015) who employ a static cooling function with a higher equilibrium value for the temperature. For the high density environments we thus observe that the dense shell pushes further into the ISM and thus increases the momentum further until it reaches pressure equilibrium with the ambient ISM. In reality, this would be prevented by the much more complex structure of the ISM (e.g. overlapping remnants or interactions with stellar wind or magnetic fields), which we will discuss in more detail in section 6. For the remnants in the lower density media the situation is different as they have an equilibrium temperature of around $10^{4}$ K in the surrounding ISM, resulting in quickly after the ST-phase established equilibrium state between shell and ambient medium, where we then fully capture the momentum conserving snow-plough phase.

In Figure 9 we show an updated version of Figure 4 of Naab & Ostriker (2017) with a comparison of the momentum gain for different ambient densities to previous estimates by Cioffi et al. (1988); Martizzi et al. (2015); Iffrig & Hennebelle (2015); Li et al. (2017); Kim & Ostriker (2015); Walch et al. (2015); Geen et al. (2016). Our results agree well in particular with Kim & Ostriker (2015) (despite a different cooling function) and show a momentum boost compared to the assumed initial SN momentum of p0 = 14181 M*⊙* km s*-1* from a factor $\sim 10$ at $n=100$ cm*-3* to a factor $\sim 50$ at $n=0.001$ cm*-3*. At low resolution this boost can be lower up to a factor of two at the highest densities. As we carried out all simulations at three different ambient metallicities, we can obtain best fit relations for all three metallicities that are given by:

[TABLE]

[TABLE]

[TABLE]

These fitting relations are shown as the red dashed-dotted ($Z=0.01Z_{\odot}$), the red solid ($Z=0.1Z_{\odot}$) and the red dashed line ($Z=1Z_{\odot}$) in Figure 9. The momentum is always measured at the end of the ST-phase. We note that a remnant can still generate momentum in the pressure driven snow-plough phase until pressure equilibrium between the bubble and the shell is reached, which can lead to higher momenta by a factor of two (Haid et al., 2016; Kim & Ostriker, 2015). Compared to other studies (e.g. Haid et al., 2016; Kim & Ostriker, 2015) we follow the momentum gain for $Z=0.1Z_{\odot}$ and $Z=0.01Z_{\odot}$. We find a weak dependence of the final terminal momenta generated by each remnant in different metallicity environments which mostly affects the terminal momenta in the high ambient density regimes (roughly a factor of 2) while the momenta in the lowest ambient environments remain unchanged.

3.5 Chemical composition

We show the chemical evolution of SN-remnants for H2, H, and H*+* in the WNM, the CNM, the CNM-n10, and the CMM in Figure 10. In the WIM (top left) the initial material is fully ionised and H*+* stays the dominating chemical component of the shell even after $2$ shell formation times where roughly $60$ per cent of the remnant are in full ionisation. However, at ten shell formation times the remnant becomes fully neutral and also forms a tiny fraction of molecular hydrogen (below a solar mass). The CNM (top right panel) initially consists of neutral hydrogen but the whole remnant is ionised by the blast wave. Only after shell formation the remnant cools and the shell is becoming 100 per cent neutral again after $10$ shell formation times. In these remnants molecular species can be ignored. For the CNM, we show a resolution study for the chemical fractions within the remnant. At the lowest mass resolutions the chemical composition of the single remnant is not captured accurately. At this resolution the remnant is not fully ionised in the beginning and initially more mass is kept in the neutral state rather than fully ionised by the blast wave as expected. The remnant is fully neutral before shell formation. The same is true for the $10$ M*⊙* runs, although we note that here we already find a slightly dominating ionisation fraction of the remnant compared to the neutral component. For the two highest resolution runs we find that the blast wave is fully ionising the ambient medium, keeping 70 to 80 per cent of the gas ionised until shell formation. After that, cooling is dominating and the remnant starts to form neutral hydrogen and reaches the fully neutral state roughly $10$ shell formation times after the appearance of the blast wave. Despite the fact that the molecular fractions are low, they also vary within at least an order of magnitude at the end of the simulation, where the low resolution remnants produce more molecular hydrogen than the high resolution ones, which leads to a more dominant cooling process through molecules at low resolution. This indicates that the non-equilibrium chemical model overestimates the cooling rate at low resolution increasing the issues with the well known cooling problem.

In higher density media (CNM-n10) we find the same behaviour (bottom left panel of Figure 10). However, we note that initially a small fraction of hydrogen of around 20 per cent is kept in the neutral phase and the formation rate of molecular hydrogen is slightly enhanced (half an order of magnitude) towards the end of the simulation. In the CMM medium (bottom left of Figure 10) we find that the some of the initially swept up mass is kept in the neutral state (again roughly 20 per cent) while the rest if fully ionised. The formation rate of molecular hydrogen increases again and at the end of the simulation we find about half an order of magnitude lower molecular fractions than ionised rest mass of the remnant. If we continued the simulation longer the built-up of molecular hydrogen would be significant. However, we assume that remnants overlap before they reach more than $30$ Sedov-times.

Finally, we briefly discuss the dynamical impact of the chemical model. We carry out one reference run where we override the evolution of the chemical abundances and only enable our non-equilibrium cooling routines. Not evolving the chemical abundances has mainly an impact on the cooling behaviour of the gas. We only find a marginal impact of the chemical model on the overall evolution of the supernova remnant on the percent level. However, the chemistry effects the initial momentum build-up and the peak hot mass in the beginning. The initial momentum is slightly increased in the run without the active chemical evolution, but converges quickly to the momentum evolution that is observed with the active chemistry. The peak hot mass is slightly increased without the chemical evolution, which can be explained by the missing cooling channel over the molecules. Furthermore, we want to point out that the radiation which is emitted by the remnant would certainly impact its chemical evolution which is a potential caveat of the model. However, the simulations are also too low in resolution to make valuable prediction for X-ray properties of supernova-remnants. To accurately investigate this, one would have to carry out higher resolution simulations with lower dimensionality like Thornton et al. (1998), Blondin et al. (1998), or more recent Badjin et al. (2016). However, We investigate a possible feedback scheme that aims for the input of the terminal momentum and the generation of the hot phase in galaxy scale simulations in a resolved fashion. Although the chemistry might be insignificant in isolated events that picture can be quickly altered if we consider multiple events that lead to overlapping higher density remnants with shorter cooling times and higher formation rates for molecular hydrogen, which will certainly influence the distribution of internal energy and therefore its dynamics within the ISM.

3.6 Comparison of different Hydro-solvers

We carried out some of our simulations employing different solvers for hydrodynamics. All simulations at our highest mass resolution are carried out with our pressure-energy SPH solver. In Figure 12 we show the comparison between the MFM-solver (solid) and the SPH-solver (dashed) for all physical quantities that we investigate. Overall, we find good agreement of the two solvers at the highest mass resolution. Most of the physical properties of the shock are captured within a few per cent. The major differences between the solvers are given in the behaviour of the thermal and kinetic energy (upper left panel of Figure 12). For the SPH-solver we find that the thermal energy is radiated away at a slightly earlier time compared to the MFM solver leading to shorter shell formation times (c.f. equation 14). Further, the SPH-solver under predicts the kinetic energy within the ST-phase. The other difference that we point out is that the SPH-solver predicts a lower temperature in the hot bubble by half an order of magnitude. Therefore, the pressure within the bubble is slightly lower, which can potentially lead to less momentum generation within the PDS phase for this method. For the SPH-solvers we evaluate the results in section 5 in more detail.

4 Build up of the hot phase

Following Naab & Ostriker (2017) we can derive the expectation value for a hot phase to be established by SN-feedback. To derive the expectation value we follow a semi-analytic approach. From our simulations we can derive the radii of the remnant as a function of the environmental density. As we do not have data for different explosion energies we lack the dependence on the explosion energy compared to similar studies (e.g. Kim & Ostriker, 2015). However, we tested different metallicities and can therefore carry out the semi-analytic calculations for different metallicity regimes. First, we need to derive the dependence of the radius on the ambient density. We can fit our results for all three different regimes. Therefore, we obtain the radius of the shell at the end of the ST-phase. We find the following relations.

[TABLE]

[TABLE]

[TABLE]

The expectation value for the hot phase of a SN-explosion that goes off within the hot phase of a previously occurred SNe can be written as follows

[TABLE]

where $S$ is the SN-rate and rst is the remnants radius at the end of the ST-phase. By substituting rst with the equations 16 to 18 and tst with equations 10 to 12 we find a strong power law dependence for the expectation value of the hot phase for all three different metallicity regimes. The expectation value for the hot phase can then be determined via:

[TABLE]

[TABLE]

[TABLE]

Thus, we see that the expectation value for the hot phase varies only very weakly with the metallicity and we can determine it for typical values within the ISM. We assume a typical surface density of the Milky Way within the solar neighbourhood around 10 M*⊙* kpc*-2*. We follow Naab & Ostriker (2017) and assume a Salpeter mass function and a disc height of 250 pc. From that we can determine the SN-rate $S$ for solar neighbourhood conditions and obtain 280 kpc*-3* Myr*-1*. By assuming the typical medium density within the Milky Way of $n=1$cm*-3* this gives a very low number for the expectation value for the hot phase in order of $10^{-5}$ to $10^{-4}$. However, for SNe in lower density environments the expectation value for the hot phase becomes quickly larger and for environmental densities below $0.01$ we find an expectation value greater than one due to the strong power law dependence of the expectation value with the density of the ambient medium. We note that this has consequences for the hot phase in the ISM. Once a configuration of the ISM is reached where SNe go off in low density environments our model predicts higher and higher expectation values with decreasing densities. SNe tend to go off in lower and lower density media. Cooling times for the gas increase alongside with the sizes of the remnants. Finally, this becomes a runaway process with dominating volume filling hot phase. This picture is in very good agreement with the findings from our simulations. Further, the model is in good agreement with the simulations carried out by Girichidis et al. (2016) who investigated the gas densities in galactic outflows and find a peak of the number density of roughly $0.01$ cm*-3*. While the remnants in the CMM have a size of a few pc, the size of remnants in the warm and hot ionised phase of the ISM can reach values up to $800$ pc. Moreover, in Figure 11 we show that the hot mass of the remnant significantly increases as a function of the environment. We measured the peak hot mass (maximum mass of the bubble) that is generated by SN-blastwaves in different environments and find strongly increasing masses in the hot phase in remnants below $0.1$ cm*-3*. We find the following best fit relations.

[TABLE]

[TABLE]

[TABLE]

Therefore, remnants in low density environments contribute much more to the build up of the hot phase by generating high expectation values for the hot phase alongside with high mass fractions of the hot gas. Remnants in high density media on the other hand have very low expectation values for a hot phase to form and generate only around $600$ M*⊙* of hot gas in the bubble which quickly cools away once the shell forms.

This has consequences for the formation of a wind driven by SNe. In this picture a wind is driven by the hot phase via the formation of a superbubble which interior is heated by SNe expanding into the bubbles of preceding SNe. Once the bubble breaks out of the disc the pressure within the bubble interior can push the shell outwards into the CGM. Hereby the outflow velocity would be limited by the sound speed within the CGM. Given a virial temperature of $10^{6}$ K for the CGM of the Milky Way this can drive outflows with a few 100 km s*-1* as they are observed (e.g. Genzel et al., 2011).

5 Effects of thermal conduction

Using our fiducial SPH-solver we investigate the behaviour of SN-remnants under the effect of thermal conduction. In supernova remnants, thermal conduction can influence the interface between the hot bubble and the shell by redistributing thermal energy from the bubble into the shell and mass from the shell into the bubble (Keller et al., 2014; El-Badry et al., 2019). With a temperature dependent thermal conduction coefficient in the Spitzer-limit (see Sec. 2.2) we investigate differences in remnant morphology and chemistry. All our highest resolution runs have been carried out including the effect of thermal conduction.

In Figure 13 we show the radial profiles of the density (top left), the temperature (top left), the velocity (bottom right) and the pressure (bottom right) for the simulation CNM-Lv4. Comparing to Figure 2, we find that the temperature in the run with thermal conduction is about an order of magnitude lower than in the case without thermal conduction, while the pressure remains roughly constant. The lower bubble temperatures result form heat flux from the hot bubble to the colder shell. The pressure is similar due to mass flux from the shell to the bubble. This effect has a general impact on the evolution of the remnants when thermal conduction is included.

Conduction describes a heat flux from the hot gas to the cold gas, which is counterbalanced by a mass flux from the cold gas to the hot gas constrained by energy conservation. This is valid before shell formation and leads to slightly higher mass fractions of the hot gas. However, in SN-remnants cooling is active after shell formation and energy conservation is violated. Therefore, the mass flux from cold to warm gas is prohibited and the energy from the bubble that enters the shell is radiated away immediately. This can impact the general physical properties of SN-remnants. Temperature is reduced and density slightly increases. While there are quantities that are only weakly affected by thermal conduction (e.g. remnants size, velocity structure), there are others which are affected significantly. First, we note that the ST-times (times of shell formation) are shorter in the presence of thermal conduction by roughly $20$ per cent. We show this in Figure 14 where we plot a comparison of the shell formation times that we find in our fiducial runs with MFM (red solid line), the SPH (red dashed line) and with SPH plus thermal heat conduction (red dahsed dotted line). Conduction leads to a redistribution of internal energy from the hot gas to the cold. This reduces temperatures and leads to a slightly enhanced cooling. Therefore, cooling times shorten and the energy conserving ST-phase terminates earlier. We find modified fitting formulas with conduction. For completeness we give the fit we obtained with SPH.

[TABLE]

[TABLE]

Because of the earlier termination of the ST-phase, the momentum input with thermal conduction is lowered compared to the runs without thermal conduction. In Figure 15 we show the momentum input as a function of the environmental density and find that the SPH solver gives lower momenta by roughly $10$ per cent while thermal conduction again reduces the terminal momentum by $20$ per cent compared to the results that we obtain with our fiducial MFM solver. We overplot the data points from the simulation as red dots (SPH) and red crosses (SPH plus thermal conduction). The modified fitting formulas for the results with SPH and conduction yield:

[TABLE]

[TABLE]

However, the largest impact can be seen in the build up of the hot mass which is reduced by 40 per cent. We show the results in 16 for MFM (solid line), SPH (dashed line) and conduction (dashed-dotted line). While the results for MFM and SPH only differ in the per cent regime, the hot mass is reduced $40$ per cent in the runs with thermal conduction. We note that before cooling takes over the hot mass slightly increased due to the mass flux from warm to hot gas. For thermal conduction we find the dependence of the hot mass as function of the environmental density as follows:

[TABLE]

Further, we find an impact of thermal conduction on the chemical composition of the SN-remnant. Overall, the effect is minor. Therefore, we show the results for the run CNM-cond in Figure 17. Initially, the majority of the gas is ionised in the case with and without conduction but with conduction the ionising fraction is lightly reduced (by roughly 5 per cent). Because cooling is slightly more efficient in the presence of conduction, the build-up of neutral hydrogen in the shell is slightly enhanced. This then leads to higher formation rates of molecular hydrogen which is increased by a factor of two.

6 Applications to astrophysical systems

In this section we briefly discuss applications to larger scale ISM and galactic scale systems to show that it is possible to use the presented feedback scheme at the target resolution of solar mass and sub-parsec resolution and achieve convergence on the momentum and hot phase generation in a global large scale simulations just by resolving the momentum and hot phase of isolated remnants. Moreover, we provide an insight into SN-remnants in structured media which can not be handled in higher resolution one dimensional SN-blast waves studies as these systems miss the radial symmetry needed for justifying a lower dimensional approach.

6.1 The supernova driven ISM

The first example we want to study is the SN-driven ISM. In Figure 18 we show the density distribution of the SN-driven ISM in a periodic box with a side length of 256 pc with solar neighbourhood conditions for the background UV-field but $0.1$ Z*⊙* which represents an ISM environment as present at high redshift. However, we adopt a lower SN-rate of 1.5 supernova per Myr and follow the evolution of the SN-driven ISM for 100 Myr. The four different panels in Figure 18 represent the four different mass resolutions after 100 Myr of evolution. Visually, we note some differences in density structure. In the lowest resolution run we observe that the resulting ISM-structure is only very badly resolved and the hot phase is mostly missing in this simulation as the canonical SN-energy falls to over cooling and is radiated away shortly after deposition without depositing momentum or generating the hot phase ans subsequently it does not provide the ISM-pressure to subsequently generate a SN-driven outflow. Finally, we investigate the SN-remnant properties that we obtain from a feedback event in a structured (SN-driven) medium. In Figure 20 we show the blast wave evolution for 4 different point of a blast wave that evolves into the SN-driven ISM for our highest resolution runs. In this scenario we generate the initial conditions by using the SN-driven ISM simulations that we show in Figure 18 and select a snapshot with a nicely build-up turbulent structure that provides an realistic environment for the explosion of massive star as we would obtain it in galactic scale simulation. The canonical SN-energy of $10^{51}$ is hereby distributed over the $32$ particles closest to the centre of the box. Thus the SN-even occurs in the centre of the box at $x=180$ pc and $y=180$ pc. The supernova goes of in the dense filamentary structure that can be observed in the upper left panel of Figure 20. With this event we represent a supernova-explosion in the CMM. As the blast wave occurs in the centre of a filament the feedback is most effective into the direction of the steepest pressure gradient. In this scenario this coincides with the steepest density-gradient. Therefore, we obtain a highly asymmetric blast wave structure with a faster expansion of the blast wave in negative and positive z-direction. While the blastwave can expand freely in negative y-direction it expands into an already bubble-like in positive y-direction and heats the material in this already present bubble even more (upper right panel of Figure 20). Despite the fact that the feedback is ineffective parallel to the filament we can still observe in the upper right panel of Figure 20 that it get disrupted, by the momentum that is deposited during the ST-phase which remains resolved in this simulation. After the end of the ST-phase (bottom panels of Figure 20) the material in the irregular formed shell is further pushed out via the pressure-support by the hot material inside the newly formed bubble, generated in the resolved ST-phase.Finally, we show the detailed blast wave properties in Figure 21 and point out the differences to the results of the runs that we obtained with a non-structured, constant density ambient medium at rest. The main difference is that in the context of an already structured ambient medium the heating by the supernova-explosion is most effective in the direction of the steepest pressure gradient. Thus the remnant becomes highly non-isotropic. While it is relatively easy to follow the momentum, the energy distribution and the hot mass fraction, of this non-isotropic remnant it is not straight forward to determine the cold mass as most of the other parts of the ISM are still cooling while they are not yet effected by the heating from the supernova. The velocity can only be determined as averaged value out of shell mass and shell momentum. We show the central quantities of the remnant in a structured medium in Figure 21 with the momentum (top left), the energy (top right), the velocity (bottom left) and the remnant mass (bottom right). We can see that the remnant, despite its irregular shaped form still undergoes the ST-phase, the PDS-phase and the MCS-phase, which can be identified in the momentum and energy distribution. Velocity and cold mass are harder to track in this context while the hot mass can be determined more easily. we note that we find a factor of 2.5 more hot mass in the structured medium compared to the remnants in homogeneous media. This can again be explained by the more efficient heating alongside the steepest density gradient.

6.2 Galaxy scale simulations

One major achievement that can be realised at the target resolution is a full galactic disc simulation that resolves the feedback of massive stars, with the focus on SN-feedback of Type II core collapse supernovae. In Figure 22 we show a simulation of an isolated dwarf galaxy that has been simulated with our code and the feedback scheme that we have tested into depth in the scope of this work to investigate the momentum and hot-phase generating properties that can be obtained with a very simple SN-feedback scheme at target resolutions where the the single feedback events’ ST-phase is resolved. This kind of simulation shows that SN-feedback can be effective once the ST-phase can be resolved in the majority of the ambient media in which the supernova occurs. With this simple, but resolved feedback scheme it becomes possible to self-consistently study ISM-turbulence and galactic scale outflows in a galaxy scale simulation, without the need for more complicated and less accurate treatments of sub-grid feedback physics, like direct momentum input or a mixed injection scheme that assumes the ST-solution by construction. For our scheme we are able to show that the vast majority of SN-events occurs in environmental densities for which our isolated blast wave studies clearly indicate convergence on momentum and hot phase generation at a mass resolution of 1 M*⊙*. Together, with the results of our presented SN-driven ISM-environments it becomes possible to also achieve hot phase-dominating Volume filling fractions in a galactic scale simulation. Thus far, this has only been achieved in small galactic patches. Finally, we note that with this kind of galaxy-scale simulation we are able to achieve one major step forward in the modelling of ISM-physics in a galactic context. Although, we treat the stars that we form as single stars by randomly drawing them from the IMF, this still represents a sub-grid model for star formation. We then follow the massive stellar population in terms of photo-ionising radiation and photo-electric heating in a Stroemgren-approximation (e.g. Hopkins et al., 2011) by mapping the density distribution with a healpix-algorithm and calculating the local column densities with treecol (Clark et al., 2012). Although, we neglect re-emission in this context, we still include an approximation for radiation of the massive stars. Once these massive stars reach the end of their life times they explode as core collapse SNe and every single remnant is able to provide its momentum and hot-phase in the ST-phase removing the sub-grid constraint for supernova-feedback from our galaxy simulations.

7 Conclusions

7.1 Summary

We carried out three dimensional simulations of isolated SNe that include a non-equilibrium cooling and a chemical model that tracks the formation (and destruction) of H2. The canonical SN-energy of $10^{51}$ erg is coupled by pure thermal injection to one kernel size the centre of a box with uniform density and tested four different mass resolutions of $0.1$ M*⊙, $1$ M⊙, $10$ M⊙* and $100$ M*⊙*. The energy is distributed weighted by the kernel. We focus on the convergence of physical properties at the end of the ST-phase to constrain the resolution requirements for a SN-feedback scheme that can be applied in high resolution simulations of galaxy formation and evolution. Specifically, we tested the behaviour of two different, widely used numerical methods for solving Euler’s equations. We tested the implications of this feedback scheme with a modified ’pressure-energy’ SPH solver and the higher order meshless finite mass (MFM) solver which utilises a second order reconstruction of Euler’s equations by solving the Riemann problem on the one dimensional surface between two fluid tracers to obtain the fluid fluxes. The usage of the chemical model allows us to follow the formation and destruction processes of molecular, neutral and ionised gas within isolated SN-remnants as they currently can be modelled in simulations of galaxy formation and evolution. We carried out reference runs with thermal conduction to investigate the effects of the heat flux between SN-bubble and the cold material in the remnant.

Morphology of the shock: We can capture three stages in the evolution of a SN-remnant. In the ST-phase we resolve the adiabatic regime defined by the Rankine-Hugoniot jump conditions with an accuracy of up to $10$ per cent that is limited by the behaviour of our underlying numerical scheme. After the ST-phase the remnants enters the PDS-phase in which they can further increase momentum due to the high pressure in the bubble behind the shell. This pushes the shell further outwards until pressure equilibrium with the shell is reached and the shell moves forward with roughly constant momentum and reaches the MCS-phase. Finally, the remnant starts to merge with the ambient ISM. We can capture the shocks morphology reasonably well at our three highest resolution levels. However, at the lowest resolution level the shock remains unresolved. 2. 2.

Dependence on the environment: We carried out explosions for six different environmental densities of SN-remnants alongside the equilibrium cooling curve. For all these remnants we computed the end of the ST-phase and find a power law scaling of the termination time of the ST-phase which is within a 30 per cent agreement compared to the work of Blondin et al. (1998), Petruk (2006) and Haid et al. (2016). We investigate eight physical quantities as a function of time including radius, momentum, energy distribution, velocity structure, remnant mass, bubble-mass, shell-mass as well as the bubble temperature. At the end of the ST-phase our resolved runs agree very well with the results of studies of the same kind (e.g. Kim & Ostriker, 2015; Haid et al., 2016), although we note that our default metallicity is $0.1$ Z*⊙* because we aim for understanding SN-feedback in low metallicity environments. Most quantities can be resolved at a resolution of $10$ M*⊙* in all relevant density regimes of the ISM. We specifically highlight this in the context of the momentum that has been generated during the ST-phase and the generation of the hot phase (the hot mass and the temperature evolution in the bubble). The feedback can then be resolved as a combination of those two quantities. Momentum is needed to move particles through the ISM and generate turbulence. The hot phase is needed to generate pressure in the ISM which is necessary to launch galactic winds that can impact the CGM of a galaxy. Further, we determine a fit to the SN-momenta that are injected in different regions of the ISM based on our highest resolution remnants which can be adopted for feedback schemes in simulations of galaxy formation and evolution that do not resolve the ST-phase explicitly. 3. 3.

Chemistry of single SN-remnants: We investigate the chemical composition in the simulated SN-remnants. We note that the remnants in the low density environments contribute to the build up of the hot phase but the remnants become fully neutral after a few cooling times. They do not contribute to the formation of molecular hydrogen because the densities in the swept up mass remain low. In the CNM the mass in the swept up region is fully ionised. After cooling sets in the remnant becomes neutral after a few cooling times again. However at the end of the simulation (t=1 Myr) we find around one per cent of the remnant mass in neutral hydrogen. In the CMM the SN destroys most of the molecular hydrogen and leads to an ionisation of the swept up mass. We find that a small fraction (around $10$ per cent) of the remnants mass remain in the neutral hydrogen phase. After cooling set in most of the material cools quickly and recombination with the free electrons leads to the formation of neutral hydrogen. 4. 4.

Role of thermal conduction: Thermal conduction leads to a decrease in temperature and to an increase in density in the shell. Thermal energy can be transported from the hot phase to the cold phase while mass from the cold phase is transported to the hot phase. Once cooling becomes relevant the mass flux towards the hot bubble is suppressed because the thermal energy that is transported from the bubble to the shell is instantly radiated away. Thermal conduction makes cooling slightly more efficient and leads to slightly less momentum generation (by $10$ per cent). The quantity that is most affected by thermal conduction is the peak hot mass within the hot bubble which decreases by roughly $40$ per cent compared to the runs without conduction. Due to an increase of the densities within the shell the formation rate of molecular hydrogen increases roughly by a factor of two. Moreover, we find good agreement of our results with the recent work of El-Badry et al. (2019) who investigated the effects of cooling and thermal conduction in simulations of one-dimensional super-bubbles. 5. 5.

Applications of the presented feedback scheme: The presented feedback scheme is perfectly suited to be used in future simulations of galaxy formation and evolution that focus on sufficient modelling of the momentum build-up and hot phase generation that can be obtained at a target resolution of a few solar masses to obtain convergence on these quantities over a dynamic range of at least six orders of magnitude. While we extensively tested the properties that can be obtained in isolated unstructured ambient media, we also investigate the case of a feedback event in a highly supersonic, turbulent medium to investigate the momentum and hot-phase build-up in a more realistic environment. Moreover, we investigate the effect of multiple SNe in a closed box model to accurately derive the Volume filling factors for the hot-phase that can be obtained with the feedback scheme. The volume filling hot phase is of importance for the explanation of large scale galactic outflows as we already discussed in section 4. We find a dominating hot and warm volume filling fraction for our highest mass resolutions while we generally find that at lower mass resolution the hot volume filling fraction remains low. Thus to self-consistently drive galactic outflows with the feedback that we present it is necessary to have a mass resolution higher than $10$ solar masses.

7.2 Model limitations

Although we find good agreement of our results with other theoretical studies we are still limited by some assumptions we already made in section 2. For the isolated remnants we assumed that they occur in an environment with a constant ambient density. Given the highly turbulent structure of the ISM (e.g. Elmegreen & Scalo, 2004) this is a rather simplified assumption. Further, we excluded the fact that massive stars normally shape their environment and ionise the surrounding ambient medium due to stellar winds prior to the core collapse event, which renders our assumption of a core collapse event in a fully molecular medium untrue. However, because in this regime we initially destroy all of the molecular hydrogen it does not influence the formation rate after shell formation. Further we neglected that each remnant ejects mass in the explosion, which initially alters the chemical composition of the remnant prior to shell formation. The injection of metals and dust in the explosion would significantly shorten the cooling time and increase the formation rate of the molecules. In a galactic context we can resolve the momentum and the hot phase in the Sedov-Taylor phase, while other properties of the a SN-remnant cannot be captured accurately. For example at the presented mass resolutions it is not possible to reach the predicted analytic value of the density in the adiabatic phase (before cooling dominates) of the remnant which can certainly change the on the fly updated cooling and molecular formation rates. Further, we do not consider magnetic fields in the shock. This could change the picture again because of the additional pressure component (the same is valid or cosmic rays). Moreover, in the case of the runs including thermal conduction it remains unclear to which degree the behaviour is driven by the numerics of the underlying scheme. Here, the crucial point is to resolve the structure of the hot and the cold phase which we do in our two highest resolution runs. If this can be done the effect is physical, if not it remains unclear whether the results have any physical impact. Further, we note that our resolution is far too low to resolve the internal heat fluxes within the bubble or the shell (once it has formed).

Acknowledgements

We thank the referee S.A.E.G. Falle for careful comments on the draft which significantly improved the quality of the paper. UPS thanks Eirini Batziou, Andreas Burkert, Chris Byrohl, Klaus Dolag, Joseph O’Leary, Rhea Silvia Remus, Felix Schulze and Simon White for helpful discussions and insights on the physics of the interstellar medium. UPS thanks Volker Springel and Prateek Sharma for useful discussion on thermal heat conduction. UPS thanks Jeong-Gyu Kim for his comments on the ISM in equilibrium. The authors gratefully acknowledge the computing time granted by the c2pap-cluster in Garching under the project number pr27mi where most of this work has been carried out. UPS and BPM are funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) with the project number MO 2979/1-1. UPS and TN acknowledge support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germanys Excellence Strategy - EXC-2094 - 390783311 from the DFG Cluster of Excellence "ORIGINS". The Flatiron Institute is supported by the Simons Foundation. SWG acknowledges support by the European Research Council via ERC Starting Grant RADFEEDBACK (no. 679852) and by the German Science Foundation via CRC956, Project C5.

Appendix A Different regions of the ISM

In this section we show the results for the runs CMM (Figure 23) and WNM (Figure 24) with the MFM-solver to investigate the solvers behaviour in capturing the ST-phase as a function of environment. For the higher density environment it is much more difficult to resolve the ST-phase at low mass resolution. While most properties are still within 50 per cent in comparison to Kim & Ostriker (2015) the hot phase remains less and less resolved at low mass resolution. The terminal momentum is still well resolved even at lower resolution. However, as we pointed out in section 3.2 this is a result of the feedback scheme where the thermal energy is injected into the $32$ nearest neighbours which leads to an overestimate of the shell momentum that scales with the particle mass. Therefore, the terminal momentum remains a weak measure to determine whether the ST-phase is resolved or not. Even, if ’correct’ terminal momentum is injected and particles move through the volume the feedback remains ineffective as long as the temperatures in shell and bubble are unresolved which then generate the pressure in the ISM. For the low density environments all important physical quantities can be properly captured with a mass resolution of 10 M*⊙*. The crucial point in this regime is the long cooling times which are longer than $0.1$ Myr. Because of this long cooling time and the negligible cooling losses it the hot phase in the bubble has time to build up and can be resolved even at lower mass resolutions. As already pointed out above the momentum seems to be converged between all mass resolutions in this environment which is again due to the fact that we sweep up too much mass with our feedback scheme within the ST-phase which counterbalances the poorly resolved velocity structure of the shell which leads to the correct momentum even at lowest mass resolution. However, because in this regime we also revolve the temperature structure of the shell and the bubble the feedback remains resolved because it generates enough pressure within the ISM.