The Global Star-Formation Law by Supernova Feedback

TL;DR

This paper presents a supernova feedback-based model explaining the origin of the Kennicutt-Schmidt relation in galactic discs, supported by simulations showing self-regulation of star formation through supernova bubbles.

Contribution

It introduces a simple, physics-based model linking supernova bubble dynamics to the star-formation law, independent of microscopic star formation recipes.

Findings

The KS relation emerges from supernova bubble self-regulation.

Simulations confirm the constant filling factor of SFR-suppressed bubbles.

Analytical and simulation results yield a KS slope around 1.5.

Abstract

We address a simple model where the Kennicutt-Schmidt (KS) relation between the macroscopic densities of star-formation rate (SFR, $\rho_{\rm sfr}$) and gas ($n$) in galactic discs emerges from self-regulation of the SFR via supernova feedback. It arises from the physics of supernova bubbles, insensitive to the microscopic SFR recipe and not explicitly dependent on gravity. The key is that the filling factor of SFR-suppressed supernova bubbles self-regulates to a constant, $f\sim 0.5$. Expressing the bubble fading radius and time in terms of $n$, the filling factor is $f \propto S\,n^{-s}$ with $s\sim 1.5$, where $S$ is the supernova rate density. A constant $f$ thus refers to $\rho_{\rm sfr} \propto n^{1.5}$, with a density-independent SFR efficiency per free-fall time $\sim 0.01$. The self-regulation to $f \sim 0.5$ and the convergence to a KS relation independent of the local SFR recipe are demonstrated in cosmological and isolated-galaxy simulations using different codes and recipes. In parallel, the spherical analysis of bubble evolution is generalized to clustered supernovae, analytically and via simulations, yielding $s \simeq 1.5 \pm 0.5$. An analysis of photo-ionized bubbles about pre-supernova stars yields a range of KS slopes but the KS relation is dominated by the supernova bubbles. Superbubble blowouts may lead to an alternative self-regulation by outflows and recycling. While the model is over-simplified, its simplicity and validity in the simulations may argue that it captures the origin of the KS relation.

Tables

Table 1. Table 1: Isolated-galaxy simulations. The local SFR recipe is ρ sfr ∝ n s loc proportional-to subscript 𝜌 sfr superscript 𝑛 subscript 𝑠 loc \rho_{\rm sfr}\propto n^{s_{\rm loc}} . The global KS slope s 𝑠 s , the rms scatter about it, and the hot filling factor f 𝑓 f are quoted for the simulations with gas fraction 15% and 50% respectively. The KS relation refers to a linear fit at n > 50 ​ cm − 3 𝑛 50 superscript cm 3 n>50\,{\rm cm}^{-3} , and f 𝑓 f is determined in the relaxed phase at t > 200 ​ Myr 𝑡 200 Myr t>200\,{\rm Myr} .

Isolated-galaxy simulations
SFR $s_{\mathrm{loc}}$	feedback	global $s$	log rms	$f$ hot
		15% 50%	15% 50%	15% 50%
standard 1.5	SN+HII	1.41 1.44	0.07 0.05	0.38 0.35
shallow   1.0	SN+HII	1.28 1.23	0.13 0.15	0.42 0.41
steep      2.0	SN+HII	1.47 1.52	0.11 0.12	0.32 0.27
standard 1.5	SN	1.34 1.39	0.06 0.06	0.33 0.31
standard 1.5	HII	1.87 1.88	0.16 0.14	0.17 0.21
standard 1.5	none	2.04 1.94	0.23 0.19	0.12 0.18

Table 2. Table 2: ν 𝜈 \nu is the number of SNe in the cluster and t b subscript 𝑡 b t_{\rm b} is the burst duration of the cluster. Δ ​ r Δ 𝑟 \Delta r is the grid cell size and r inj subscript 𝑟 inj r_{\rm inj} is the radius of the region where the SN energy is injected.

Spherical simulations of SN explosions
Case	$\nu$	$t_b$	$\Delta r$	$r_{\mathrm{inj}}$
		$\mathrm{Myr}$	$\mathrm{pc}$	$\mathrm{pc}$
single SN	1	0	0.024-0.0012	0.5
A1	10	0.001	0.244-0.0122	2.0
A2	10	$0.5t_c$	0.244-0.0244	2.0
B1	100	1	0.244-0.0244	2.0
C	100	$\left(33-75\right)t_f$	0.122	2.0
cloud of clusters	10, 66	30, 1	0.244-0.122	2.0

Table 3. Table 3: Global properties of the VELA 3 galaxies. The quantities are quoted at z = 2 𝑧 2 z=2 . M v subscript 𝑀 v M_{\rm v} is the total virial mass. The following four quantities are measured within 0.2 ​ R v 0.2 subscript 𝑅 v 0.2R_{\rm v} : M s subscript 𝑀 s M_{\rm s} is the stellar mass, M g subscript 𝑀 g M_{\rm g} is the gas mass, SFR is the star formation rate, and R e subscript 𝑅 e R_{\rm e} is the half-stellar-mass radius. The disc outer volume, as defined in Mandelker et al . ( 2014 ) , is given by R d subscript 𝑅 d R_{\rm d} and H d subscript 𝐻 d H_{\rm d} , the disc radius and half height. The following three quantities refer to the thin disc as analyzed here, a cylinder of radius 0.8 ​ R d 0.8 subscript 𝑅 d 0.8R_{\rm d} and half-height 0.25 ​ kpc 0.25 kpc 0.25\,{\rm kpc} : n thin subscript 𝑛 thin n_{\rm thin} is the mean density of cold gas ( T < 3 × 10 4 𝑇 3 superscript 10 4 T<3\times 10^{4} K), SFR thin is the star formation rate, and f hot , thin subscript 𝑓 hot thin f_{\rm hot,thin} is the volume filling factor of hot gas ( T > 3 × 10 4 𝑇 3 superscript 10 4 T>3\times 10^{4} K). The simulations of most galaxies end near z ∼ 1 similar-to 𝑧 1 z\sim 1 , but galaxies V16-19,31,and 35 do not reach z = 2 𝑧 2 z=2 . a in , disc subscript 𝑎 in disc a_{\rm in,disc} and a fin , disc subscript 𝑎 fin disc a_{\rm fin,disc} are the earliest and latest cosmological expansion factors within which the galaxy has a cold disc with R d / H d > 4 subscript 𝑅 d subscript 𝐻 d 4 R_{\rm d}/H_{\rm d}>4 . Galaxies V02,03,05,18,34, and 35 do not have a disc phase within the simulated period.

Properties of the VELA galaxies
Galaxy	$M_v$	$M_s$	$M_g$	SFR	$R_e$	$R_d$	$H_d$		$n_{\mathrm{thin}}$	SFRthin	$f_{\mathrm{hot},\mathrm{thin}}$		$a_{\mathrm{in},\mathrm{disc}}$	$a_{\mathrm{fin},\mathrm{disc}}$
	${10}^{12}{M}_{\odot }$	${10}^{10}{M}_{\odot }$	${10}^{10}{M}_{\odot }$	$M_{\odot }/\mathrm{yr}$	kpc	kpc	kpc		${\mathrm{cm}}^{-3}$	$M_{\odot }/\mathrm{yr}$
V01	0.16	0.22	0.24	2.65	0.93	5.15	2.57		0.115	1.39	0.69		0.38	0.38
V02	0.13	0.19	0.31	1.83	1.81	6.37	3.57		0.039	0.26	0.74		-	-
V03	0.14	0.43	0.18	3.74	1.41	5.21	2.34		0.212	2.92	0.58		-	-
V04	0.12	0.10	0.17	0.48	1.73	5.71	2.79		0.036	0.12	0.49		0.48	0.48
V05	0.07	0.10	0.15	0.59	1.81	5.36	1.98		0.014	0.11	0.60		-	-
V06	0.55	2.22	0.52	20.63	1.05	2.53	0.42		4.155	16.33	0.61		0 .17	0.33
V07	0.90	6.37	1.98	26.66	2.85	12.59	2.06		0.143	4.67	0.65		0.20	0.54
V08	0.28	0.36	0.32	5.76	0.74	4.03	1.53		0.542	4.58	0.81		0.45	0.57
V09	0.27	1.07	0.49	3.94	1.74	7.34	2.12		0.121	0.91	0.57		0.29	0.40
V10	0.13	0.64	0.21	3.27	0.46	4.51	1.19		0.396	2.07	0.59		0.27	0.56
V11	0.27	1.02	0.86	17.18	2.14	8.34	5.08		0.115	1.85	0.79		0.29	0.46
V12	0.27	2.06	0.30	2.90	1.13	6.53	1.72		0.079	0.65	0.61		0.20	0.44
V13	0.31	0.96	1.34	21.20	2.48	9.74	4.75		0.088	1.47	0.66		0.36	0.39
V14	0.36	1.40	0.79	27.50	0.32	-	-		-	-	-		0.37	0.41
V15	0.12	0.56	0.23	1.71	1.07	6.26	1.08		0.096	0.82	0.57		0.30	0.51
V16	-	-	-	-	-	-	-		-	-	-		0.14	0.24
V17	-	-	-	-	-	-	-		-	-	-		0.15	0.31
V18	-	-	-	-	-	-	-		-	-	-		-	-
V19	-	-	-	-	-	-	-		-	-	-		0.14	0.29
V20	0.53	3.92	0.73	7.26	1.72	9.57	2.75		0.053	1.09	0.76		0.11	0.44
V21	0.62	4.28	0.92	9.80	1.73	9.48	1.18		0.191	2.78	0.56		0.20	0.49
V22	0.49	4.57	0.26	12.08	1.31	4.70	0.40		1.270	9.41	0.58		0.16	0.50
V23	0.15	0.84	0.32	3.37	1.16	6.28	1.54		0.167	2.03	0.60		0.33	0.46
V24	0.28	0.95	0.49	4.39	1.68	7.29	1.95		0.071	0.82	0.60		0.37	0.48
V25	0.22	0.76	0.17	2.31	0.73	5.70	0.82		0.173	2.09	0.62		0.32	0.50
V26	0.36	1.63	0.44	9.66	0.74	5.42	1.30		0.472	4.39	0.69		0.26	0.50
V27	0.33	0.90	0.90	8.69	1.98	9.16	4.97		0.085	3.57	0.74		0.37	0.50
V28	0.20	0.27	0.38	5.72	2.32	5.66	2.97		0.160	2.74	0.74		0.41	0.50
V29	0.52	2.67	0.54	18.75	1.89	7.46	0.97		0.270	6.37	0.78		0.26	0.50
V30	0.31	1.71	0.66	3.85	1.43	9.32	1.67		0.051	0.72	0.61		0.19	0.33
V31	-	-	-	-	-	-	-		-	-	-		0.17	0.19
V32	0.59	2.74	0.60	15.00	2.58	4.98	1.06		0.822	4.27	0.60		0.16	0.33
V33	0.83	5.17	0.59	32.74	1.23	4.59	0.88		1.384	17.02	0.72		0.19	0.39
V34	0.52	1.73	0.67	14.69	1.84	5.29	1.87		0.629	6.03	0.68		-	-
V35	-	-	-	-	-	-	-		-	-	-		-	-

Equations

$${\rho }_{\mathrm{sfr}}={ϵ}_{\mathrm{ff}}\frac{{\rho }_{\mathrm{g}}}{t_{\mathrm{ff}}},$$

$${\rho }_{\mathrm{sfr}}\simeq 0.66\times 10^{-2}{M}_{\odot }{\mathrm{yr}}^{-1}{\mathrm{kpc}}^{-3}{ϵ}_{-2}{n}_0^{3/2}.$$

$${\rho }_{\mathrm{sfr}}\propto n^s,$$

$$f\propto S{n}^{-s},$$

$${\rho }_{\mathrm{sfr}}\propto S\propto n^s.$$

$$R_{\mathrm{s}}=1.15{\left(\frac{E{t}^2}{\rho }\right)}^{1/5}=5.1\mathrm{pc}{e}_{51}^{1/5}{n}_0^{-1/5}{t}_3^{2/5},$$

$$V_{\mathrm{s}}=\frac{2}{5}\frac{R_{\mathrm{s}}}{t}=1950\mathrm{km}{\mathrm{s}}^{-1}{e}_{51}^{1/5}{n}_0^{-1/5}{t}_3^{-3/5},$$

$$T_{\mathrm{s}}=\frac{3}{16}\frac{\mu m_{\mathrm{p}}}{k_{\mathrm{B}}}{V}_{\mathrm{s}}^2=5.3\times 10^7\mathrm{K}{e}_{51}^{2/5}{n}_0^{-2/5}{t}_3^{-6/5},$$

$$\Lambda \simeq \lambda (Z)T_6^{-0.7},\lambda (Z_{\odot })\simeq 1.1\times 10^{-22}\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^3.$$

$$\dot{E}(t)=-{\int }_0^{R_{\mathrm{s}}(t)}\Lambda [T(r^{\prime })]n_e(r^{\prime })n_H(r^{\prime })4\pi r^{\prime 2}d{r}^{\prime },$$

$$\Delta E(t)={\int }_0^t\dot{E}(t^{\prime })d{t}^{\prime }.$$

$$t_{\mathrm{c}}\simeq 4.93\times 10^4\mathrm{yr}{e}_{51}^{0.22}{n}_0^{-0.55},$$

$$R_{\mathrm{c}}\simeq 24.4\mathrm{pc}{e}_{51}^{0.29}{n}_0^{-0.42}.$$

$$V_{\mathrm{c}}=188\mathrm{km}{\mathrm{s}}^{-1}{e}_{51}^{0.066}{n}_0^{0.13},$$

$$T_{\mathrm{c}}=4.86\times 10^5\mathrm{K}{e}_{51}^{0.13}{n}_0^{0.26}.$$

$$R_{\mathrm{s}}\simeq R_{\mathrm{s}}(t_{\mathrm{c}}){\left(\frac{t}{t_{\mathrm{c}}}\right)}^{2/7},$$

$$V_{\mathrm{s}}\simeq \frac{2}{7}\frac{R_{\mathrm{s}}}{t}.$$

$$t_{\mathrm{f}}\simeq 1.87\mathrm{Myr}{e}_{51}^{0.32}{c}_1^{-7/5}{n}_0^{-0.37},$$

$$R_{\mathrm{f}}\simeq 69.0\mathrm{pc}{e}_{51}^{0.32}{c}_1^{-2/5}{n}_0^{-0.37}.$$

$$f_0=V_{\mathrm{sn}}{n}_{\mathrm{sn}}=\frac{4\pi }{3}{R}_{\mathrm{f}}^{3}S{t}_{\mathrm{f}},$$

$$f_0=0.24e_{51}^{1.26}{c}_1^{-2.6}{S}_{-4}{n}_0^{-1.48}.$$

$${\rho }_{\mathrm{sfr}}\propto S\propto n^{1.48},$$

$${\rho }_{\mathrm{sfr}}={\mu }_{\mathrm{sn}}S,$$

$${\rho }_{\mathrm{sfr}}=1.26\times 10^{-2}{M}_{\odot }{\mathrm{yr}}^{-1}{\mathrm{kpc}}^{-3}{\mu }_2{e}_{51}^{-1.26}{c}_1^{2.6}{f}_0{n}_0^{1.48}.$$

$${\rho }_{\mathrm{sfr}}\simeq 0.66\times 10^{-2}{M}_{\odot }{\mathrm{yr}}^{-1}{\mathrm{kpc}}^{-3}{n}_0^{1.48}.$$

$$f=({\varphi }^2+f_0^2{)}^{-1/2}{f}_0,\varphi \simeq 1.375.$$

$${\rho }_{\mathrm{sfr}}\propto {\rho }_{\mathrm{g}}^{s_{\mathrm{loc}}}.$$

$$L=\frac{\gamma E}{t_{\mathrm{s}}}.$$

$${\ell }_{38}=0.32e_{51}{t_{\mathrm{s}}}_6^{-1},$$

$$R_{\mathrm{s}}=0.88{\left(\frac{L{t}^3}{\rho }\right)}^{1/5}=76.55\mathrm{pc}{\ell }_{38}^{1/5}{n}_0^{-1/5}{t}_6^{3/5},$$

Full text

††pagerange: The Global Star-Formation Law by Supernova Feedback–A.2††pubyear: 2002

1 Introduction

The global  Kennicutt-Schmidt relation commonly refers to the observed correlation between the surface densities of star formation rate ($\Sigma_{\rm sfr}$) and molecular gas ($\Sigma_{\rm g}$), either in galactic discs as a whole or in macroscopic regions within the discs (Kennicutt, 1998; Daddi et al., 2010b). The quantities are averaged on scales larger than the disc height, namely from $\sim 100\,{\rm pc}$ to several kiloparsecs, where the average gas number densities are $n\sim 1\,{\rm cm}^{-3}$. The global relation is typically $\Sigma_{\rm sfr}\propto\Sigma_{\rm g}^{s}$, with the slope ranging about $s\simeq 1.5$ in the range $1-2$, as summarized in §8.5. This global relation may or may not be related to the local, microscopic relation between the densities on the scales of the star-forming regions, typically smaller than $10\,{\rm pc}$, where the number densities are $n\sim 10^{2-4}\,{\rm cm}^{-3}$.

Different galaxy types and sub-galactic regions, at different redshifts and environments and on different scales, may show somewhat different KS relations, which makes the overall compiled relation look rather loose. However, it has been demonstrated (Krumholz, Dekel & McKee, 2012) that the local correlation becomes particularly tight and universal once $\Sigma_{\rm g}$ is replaced by $\Sigma_{\rm g}/t_{\rm ff}$, where $t_{\rm ff}$ is the proper free-fall time in the relevant star-forming regions. They argued that this is consistent with a universal local 3D star-formation law,

[TABLE]

where $\rho_{\rm sfr}$ is the star-formation rate (SFR) density and $\rho_{\rm g}$ is the local molecular gas density, averaged over the star-forming molecular cloud.

The microscopic SFR efficiency, $\epsilon_{\rm ff}$, appears to be constant, independent of density, at a value on the order of $\epsilon_{\rm ff}\sim 0.01$. One line of evidence for this comes from direct measurements in individual resolved clouds (Krumholz, Dekel & McKee, 2012; Salim, Federrath & Kewley, 2015; Heyer et al., 2016; Vutisalchavakul, Evans & Heyer, 2016; Leroy et al., 2017). The second line of evidence comes from the correlation of SFR with HCN luminosity, which constrains $\epsilon_{\rm ff}$ because HCN emission comes from gas at a known density (Krumholz & Tan, 2007; García-Burillo et al., 2012; Usero et al., 2015; Onus, Krumholz & Federrath, 2018).111One should note, however, that there are conflicting indications for variations in $\epsilon_{\rm ff}$ (Murray, 2011; Lee, Miville-Deschênes & Murray, 2016), which are argued to be due to a bias in the methodology (Leroy et al., 2017; Krumholz et al., 2017).

To make eq. (1) more quantitative, the free-fall time can be expressed in terms of the density as $t_{\rm ff}=[32\,G\rho/(3\pi)]^{-1/2}$, where we assume that the density is dominated by the gas. With a nucleon number density $n=n_{H}=\rho/(\mu_{H}m_{\rm p})$ (we adopt $\mu_{H}=1.355$ for Solar metallicity), and denoting $n=1\,\,{\rm cm}^{-3}\,n_{0}$ and $\epsilon_{\rm ff}=0.01\epsilon_{-2}$, eq. (1) becomes

[TABLE]

It seems that this kind of relation can be extrapolated from the small scales of the star-forming clouds to larger macroscopic scales, and it can be considered as the macroscopic KS relation in 3D, the origin of which we wish to understand. In particular, we wish to figure out the origin of the slope $s\simeq 1.5$ in the macroscopic relation

[TABLE]

which is equivalent to asking why the global $\epsilon_{\rm ff}$ is not varying as a function of $n$. We also wish to figure out the origin of the amplitude of the macroscopic relation, namely what determines the global value of $\epsilon_{\rm ff}\sim 0.01$, which seems to resemble the local value.

There are in general two types of theoretical attempts to understand the KS relation (see Krumholz, 2017). The bottom-up approach, applied to the local SFR law, is based on the expected probability distribution function (PDF) of gas densities in molecular clouds and the assumption that stars form only above a constant threshold density (e.g. Padoan et al., 2014; Burkhart, 2018). The characteristic density PDF in supersonic turbulence, largely based on simulations, is expected to be a lognormal distribution, while self-gravity is expected to generate a power-law high-density tail. Integrating above the threshold density yields the microscopic $\epsilon_{\rm ff}$, which may be subject to an uncertainty in the value and constancy of the threshold density.

The top-down approach, applied macroscopically, attempts to model the KS relation as a result of self-regulation by the combined effects of gravity, accretion, star formation and feedback. These models may attempt to study the evolution of molecular clouds in a realistic inter-stellar medium. They commonly consider the balance between the momentum provided by feedback and the vertical self-gravity. They address gravitational disk instability and the driving of turbulence by feedback, by internal inflow in the disc driven by disk instability, and by accretion into the disk (e.g. Dekel, Sari & Ceverino, 2009; Ostriker & Shetty, 2011; Faucher-Giguère, Quataert & Hopkins, 2013; Krumholz et al., 2017). These models typically assume a form for the SFR law, e.g., eq. (1) with $\epsilon_{\rm ff}$ a constant independent of density, but they usually do not attempt to explain why it is so.

Here we explore the possibility that the macroscopic KS relation is naturally driven by self-regulation of the SFR by supernova (SN) feedback, insensitive to the specific small-scale star-formation recipe, and not even involving gravity in an explicit way. The simple key hypothesis is that the mass filling factor of the gas in which the SFR is suppressed by feedback self-regulates to a constant value on the order of $f\sim 0.5$. The SFR is suppressed (boosted) when $f$ exceeds (falls short of) this attractor value. We assume that a proxy for this filling factor is the volume filling factor of hot gas in supernova bubbles when they fade away. If the final bubble radius and fading time can be expressed as power laws of the ISM gas density $n$, then the hot volume filling factor can be expressed as

[TABLE]

where $S$ is the SN rate density (e.g. §2.2). In this case a constant $f$ would automatically imply a macroscopic star-formation law of the desired form in eq. (3),

[TABLE]

Indeed, as we will see in §2 based on the standard evolution of isolated single-SN bubbles (e.g., Draine, 2011, chapter 39), the predicted slope is $s=1.48$, suspiciously close to the macroscopic star-formation law, eq. (3). Seed ideas along similar lines have been proposed in earlier work (Dekel & Silk, 1986; McKee & Ostriker, 1977; Silk, 1997).

The weak dependence of the global KS relation on the assumed local SFR recipe was indicated in hydro simulations (Hopkins, Quataert & Murray, 2011; Hopkins, Narayanan & Murray, 2013). Here we explore how this self-regulation is materialized through the hypothesis that the hot volume filling factor of SN bubbles is self-regulated into a constant value. This key hypothesis will be tested below using isolated and cosmological simulations of galaxies, confirming the insensitivity of the global relation to the local SFR recipe and the dominant role of SN feedback in it. In parallel, we will analytically compute the hot filling factor as a function of ISM density for different sequences of co-local SNe in star-forming clusters, which for a constant hot filling factor will provide predictions for the KS relation. These analytic results will be tested and refined using simple spherical simulations.

Our analytic modeling makes several simplifying assumptions, including the following (to be discussed below, especially in §7 and §8):

$\bullet$ The medium outside the SN bubbles is uniform, ignoring the complexities associated with the multi-phased ISM and the origin of molecular hydrogen for star formation.

$\bullet$ Supernova feedback is negative, such that SFR is suppressed in the gas that has been swept by the SN bubbles (discussed in §8.1).

$\bullet$ The bubbles are largely confined to the galactic disc, while the possible effects of super-bubble blow-out are discussed in §8.4.

$\bullet$ The SN bubbles overwhelm the photo-ionized bubbles about the pre-SN O/B stars. This is argued analytically in §7, and demonstrated in simulations in §4.

The paper is organized as follows. In §2 we address the idealized case of randomly distributed single SNe, where in §2.1 we summarize the standard evolution of a single SN, and in §2.2 we introduce the concept of self-regulated hot filling factor, compute it for single SN bubbles and derive the KS relation. In §3 we test the validity of self-regulation into a constant hot filling factor using ART hydro-gravitational simulations of discs in a cosmological setting. The simulations are elaborated on in §A. In §4, using RAMSES simulations of isolated galaxies, we reproduce the self-regulation to a constant filling factor and the generation of a global KS relation, and show that it is insensitive to the local SFR recipe and is determined by SN feedback. In §5 we return to analytic modeling, derive the evolution of a co-local multiple SNe in different time sequences, and obtain the associated filling factor and KS relation. In §6 we use spherical simulations to test and modify the analytic predictions for such clustered SNe. In §7 we address the alternative of photo-ionized bubbles. In §8 we discuss our modeling, where in §8.1 we address the assumption of negative feedback, in §8.2 we comment on the dominance of molecular hydrogen, in §8.3 we refer to the relevance of self gravity, and in §8.4 we comment on the effects of super-bubble blowout. In §9 we summarize our results and discuss our conclusions.

2 KS Relation - Isolated Supernovae

2.1 Standard SN-Bubble Evolution

We first summarize the standard evolution of single spherical SN bubbles (e.g., Draine, 2011, chapter 39). One assumes that a supernova of energy $E=10^{51}\,{\rm erg}\,e_{51}$ explodes in a uniform medium of Hydrogen number density $n$ and with a speed of sound (or turbulence velocity dispersion) $c=10\,{\rm km}\,{\rm s}^{-1}c_{1}$.

2.1.1 The Sedov-Taylor Phase and Cooling

After a free expansion phase, dominated by the mass of the ejecta, the SN bubble enters the Sedov-Taylor adiabatic phase, where it is approximated as a point explosion ejecting energy into a cold medium of uniform density, neglecting radiative losses, the mass of the ejecta and the pressure in the medium. Based on dimensional analysis, the shock radius, velocity and temperature are

[TABLE]

[TABLE]

[TABLE]

where the time is $t=10^{3}\,{\rm yr}\,t_{3}$. The internal profiles within the bubble are assumed to obey the Sedov-Taylor similarity solution and are obtained numerically.

Cooling just behind the shock front eventually makes the bubble leave the adiabatic phase and enter the radiative phase. In order to estimate the cooling time $t_{\rm c}$, the cooling function in the relevant temperature range is idealized by

[TABLE]

This is a fair approximation for solar-metallicity gas, our fiducial case below, at temperatures in the range $10^{5}-10^{7.3}$K. The cooling rate is obtained by spatial integration over the bubble,

[TABLE]

where $n_{e}=\rho/(\mu_{e}m_{\rm p})$ and $n_{H}=\rho/(\mu_{H}m_{\rm p})$ (with $\mu_{e}=1.15$ and $\mu_{H}=1.355$ for solar metallicity). The integral over $T^{-0.7}(r)\,n^{2}(r)\,r^{2}$ is evaluated numerically for the Sedov-Taylor similarity solution. The energy loss by time $t$ is

[TABLE]

The SN bubble ends its adiabatic phase, e.g., having lost one third of its energy to radiation, after a time $t_{\rm c}$ when the shell is at a radius $R_{\rm c}$ given by

[TABLE]

[TABLE]

The shock velocity and temperature are then

[TABLE]

[TABLE]

2.1.2 The Snow-plow Phase and Fading

At $t>t_{\rm c}$, after a significant fraction of the original SN energy has been lost to radiation, a dense shell of cold gas is pushed by the pressure of the enclosed hot central volume. The mass of the dense shell increases as it sweeps up the ambient gas, and it slows down accordingly,

[TABLE]

[TABLE]

At the beginning of the snow-plow phase $V_{\rm s}\sim 150\,{\rm km}\,{\rm s}^{-1}$ propagating into a medium of T\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}10^{4}K, namely it is a strong shock.

The snow-plow phase fades away when $V_{\rm s}$ becomes comparable to the speed of sound of the ambient medium $c$. This occurs after a fading time $t_{\rm f}$, leaving behind a bubble of fade-away radius $R_{\rm f}$,

[TABLE]

[TABLE]

The mass affected by the SN bubble, which is mostly in the broadened, faded away shell, is the mass that was initially within the volume encompassed by $R_{\rm f}$. As discussed in §1, we assume that the SFR in the shell is suppressed by mechanical effects.

2.1.3 Weak Dependence on Metallicity

Slower cooling because of a lower metallicity $Z$ would make the final bubble larger. For a lower metallicity, the factor $\lambda(Z)$ in eq. (9) is smaller. For example, for $Z/Z_{\odot}\simeq 0.1$ the value drops by a factor of $\simeq 7$ (e.g. Draine, 2011, Fig. 34.1). In the expressions above, one obtains $t_{\rm c}\propto\lambda^{-1/3}$ and $R_{\rm c}\propto\lambda^{-2/15}$. This leads to $R_{\rm f}\propto t_{\rm f}\propto\lambda^{-0.053}$. The volume of a bubble at fading is then $V_{\rm sn}\propto\lambda^{-0.21}$. For example, if $Z=0.1\,Z_{\odot}$, for which $\lambda\simeq 0.14\,\lambda(Z_{\odot})$, the volume per SN bubble becomes larger by about 50%. This is a relatively weak dependence on metallicity, which we will not deal with here.

Small changes in the power-law fit to the $T$ dependence of the cooling curve, eq. (9), lead to small changes in the power indices of $n$. For example, Dekel & Silk (1986) assumed $\Lambda\propto T^{-1}$ (instead of $T^{-0.7}$), relevant for a gas with lower metallicity, and obtained at the end of the radiative phase a somewhat larger bubble with $t_{\rm c}\propto n^{-0.53}$ and $R_{\rm c}\propto n^{-0.41}$ compared to the powers of $-0.55$ and $-0.42$ in eq. (12) and eq. (13). These are again relatively small effects, which we ignore here.

2.2 Bubble Filling Factor and the KS Relation

We now derive the KS relation from the bubble fading radius and time and their dependence on $n$.

2.2.1 Constant filling factor as an attractor

Consider first the idealized case where the supernovae occur at random positions within a uniform galactic disc (of a constant height, say), with a SN rate per unit volume $S=10^{-4}\,{\rm kpc}^{-3}\,{\rm yr}^{-1}\,S_{-4}$.

The volume filling factor of hot bubbles of radius $R_{\rm f}$ after a fade-away time $t_{\rm f}$, tentatively neglecting overlaps between bubbles, is

[TABLE]

where $V_{\rm sn}$ is the volume of each bubble and $n_{\rm sn}$ is the number density of pre-fading bubble centers. Substituting $R_{\rm f}$ and $t_{\rm f}$ from eq. (19) and eq. (18) one obtains

[TABLE]

Thus, a filling factor of a fixed value of order one half defines a critical line in the $S\!-\!n$ diagram, $S\propto n^{1.48}$. This has been demonstrated to be approximately reproduced in simulations (Li et al., 2015). For a system above that line, most of the ISM is “hot” within the bubbles, where SFR is suppressed, while below the line it is unperturbed, cold and available for star formation.

Our main point here is that the natural proportionality of the SN rate and the SFR, for a given stellar initial mass function (IMF), makes this line an attractor. When above the line, most of the ISM is “hot”, the SFR is suppressed, so $S$ is suppressed and decreases down toward the critical line. When below the line, most of the ISM is cold, the SFR is free to grow, so $S$ increases toward the critical line. The galaxies are thus expected to populate a line

[TABLE]

reminiscent of the KS law.

2.2.2 The KS Relation for isolated SNe

In order to obtain the normalization of the KS relation one should translate S to SFR density, given the IMF,

[TABLE]

where $\mu_{\rm sn}=100\,\mu_{2}\,M_{\odot}$ is the mass in forming stars that generate one supernova, namely the ratio of SFR to the SN rate. We obtain from eq. (21)

[TABLE]

Using the Milky Way as an example, the speed of sound is $c\simeq 6\,{\rm km}\,{\rm s}^{-1}$, the global SFR is ${\rm SFR}\simeq 2M_{\odot}\,{\rm yr}^{-1}$ (Chomiuk & Povich, 2011; Licquia & Newman, 2015) and the SN rate is $1/60\,{\rm yr}^{-1}$, yielding $\mu_{\rm sn}\simeq 120M_{\odot}$. With $e_{51}=1$ and $f_{0}=0.5$ we obtain (for isolated SNe, ignoring overlaps)

[TABLE]

Eq. (25) coincides with eq. (2) representing the observed KS relation for $\epsilon_{\rm ff}\simeq 0.016$. This simple analysis predicts a universal relation of the form in eq. (1) with a constant $\epsilon_{\rm ff}$ independent of $n$ and a value in the observed ball park, though with a non-negligible dependence on the sound speed in the ISM.

2.2.3 Estimated Correction for Overlaps

The filling factor as computed above, $f_{0}=n_{\rm sn}V_{\rm sn}$, ignored overlaps between bubbles. If bubbles overlap, the actual volume filling factor $f$ is smaller. At low filling factors $f$ approximately coincides with $f_{0}$, but at large filling factors $f_{0}$ may become a severe overestimate of the actual $f$. In order to estimate $f(f_{0})$, we tentatively consider a random distribution of bubble centers in the volume. Figure 1 shows the result of a numerical experiment, which is well fitted by the function

[TABLE]

This gives for example $f_{0}\simeq 0.60,0.79,1.03$ for $f=0.4,0.5,0.6$ respectively. While $f$ is limited from above by unity, $f_{0}$ could in principle be larger than unity. With $f\sim 0.5$, the effect of overlap is on the order of 50%, keeping the filling factor in the same ball park. Eq. (26) can be used to correct $f_{0}$ into $f$, especially for larger filling factors. Recall, however, that this is a crude approximation for single SNe at random positions. In reality, the SNe are clustered (§5), and the correction for overlap between super-bubbles should be recalculated accordingly.

3 Cosmological Simulations: a Constant Hot Filling Factor

Before generalizing the analytic estimates to cases of clustered SNe and testing the analytic models with spherical simulations, we turn to galactic discs in full hydro-gravitational simulations that incorporate star formation and supernova feedback, both in cosmological simulations (this section) and in isolated galaxies (next section). Our main goal here is to explore the validity of the key hypothesis of self-regulation by feedback into a constant hot volume filling factor. Using the isolated-galaxy simulations we will also explore the robustness to different local recipes for star formation, and the dominance of SN feedback.

We first utilize the suite of 35 VELA zoom-in cosmological simulations. Its relevant characteristics are mentioned here, while more details are provided in §A and in references therein. The simulations are based on an Adaptive Refinement Tree (ART) code (Kravtsov, Klypin & Khokhlov, 1997; Ceverino & Klypin, 2009). The suite consists of 35 galaxies that were evolved to $z\sim 1$, with a unique maximum spatial resolution ranging from $17.5$ to $35\,{\rm pc}$ at all times. The dark-matter halo masses range from $10^{11}$ to $10^{12}M_{\odot}$ at $z=2$.

The local SFR recipe allows stochastic star formation in grid cells where the gas temperature is below $10^{4}$K and the gas density is above a threshold of $1\,{\rm cm}^{-3}$. If we attempt to approximate the local stochastic star formation recipe by an expression of the sort $\rho_{\rm sfr}=\epsilon_{\rm ff}\rho_{\rm g}/t_{\rm ff}$ (eq. 1), the efficiency would be $\epsilon_{\rm ff}\sim 0.02$. Being close to the desired global KS relation by construction through the local SFR recipe, these simulations by themselves do not explore the predicted insensitivity of the KS relation to the local SFR recipe.

Supernova feedback is implemented as a local injection of thermal energy (Ceverino & Klypin, 2009; Ceverino, Dekel & Bournaud, 2010; Ceverino et al., 2012). The energy from SN explosions (and stellar winds) is released at a constant heating rate over the $40\,{\rm Myr}$ following the formation of the stellar particle, comparable to the age of the least massive star that explodes as a type-II, core collapse supernova. A velocity kick of $\sim 10\,{\rm km}\,{\rm s}^{-1}$ is applied to $30\%$ of the newly formed stars to mimic the effect of runaway from the densest region where the cooling is rapid. The later effects of type-Ia supernovae are also included. Naturally, the $\sim 25\,{\rm pc}$ grid does not resolve the main phases of the SN-bubble evolution, making the treatment of SN feedback rather approximate.

In addition, radiation-pressure from massive stars is implemented at a moderate level with no infrared trapping (Ceverino et al., 2014). This is incorporated through the addition of a non-thermal pressure term in cells neighboring massive star particles younger than $5\,{\rm Myr}$ and whose column densities exceed $10^{21}\,{\rm cm}^{-2}$.

Galactic discs are selected for analysis from all the snapshots available in the redshift range $z=5.6-1$ in time intervals of $\sim 100\,{\rm Myr}$. The selection criterion for a disc is that the cold-gas ($T<3\times 10^{4}$K) axial ratio is $R_{\rm d}/H_{\rm d}>4$, where $R_{\rm d}$ and $H_{\rm d}$ are the disc radius and half-height as defined in Mandelker et al. (2014), see §A, yielding 25 galaxies with long periods as discs. The radii and half-height $R_{\rm d}$ and $H_{\rm d}$ at $z=2$ span the ranges $2.5-12.6\,{\rm kpc}$ and $0.4-2.1\,{\rm kpc}$ respectively, but these rather thick cylinders may include regions off the main bodies of the discs, given that many discs are warped or asymmetric. The analysis here is conservatively confined to central thin cylinders of radii $0.8R_{\rm d}$ and height $\pm 0.25\,{\rm kpc}$ as representing the main bodies of the gas discs (Fig. 2). When macroscopic sub-volumes are desired, each disc is divided into nine patches as in Fig. 2, consisting of eight equal orthogonal patches covering the ring between $0.27R_{\rm d}$ and $0.8R_{\rm d}$ and one circular central patch of radius $0.27R_{\rm d}$ such that it has the same area as the other patches.

Our main result from the cosmological simulations is presented in Fig. 3, which shows the evolution of the hot volume filling factor in the whole disc in all snapshots of all galaxies in their discy phase, where “hot” refers to $T>3\times 10^{4}$K. The snapshots of each galaxies are connected by a line of a random color. We learn that the filling factor for each galaxy oscillates about a self-regulated fixed value, roughly $f\simeq 0.6\pm 0.07$, and that this is similar for all galaxies. We find that the filling factor is largely independent of redshift and of galaxy mass. There may be an apparent weak decline of $f$ with time, and an apparent slight increase of $f$ with the mass ranking at a fixed redshift (not explicitly shown in this figure), but the significance of these trends are questionable. This supports our basic ansatz that the feedback self-regulates the bubble volume filling factor in the discs into a roughly constant value of order one half.

The choice of $T>3\times 10^{4}$K in Fig. 3 was rather arbitrary. To test for robustness, Fig. 4 shows the cumulative volume-weighted distribution of $T$ in the disc, namely the volume filling factor for gas $>T$ as a function of $T$, stacked for all the discy snapshots in the redshift range $z=5.6-1.0$. Also shown (dashed red) is the same but for the high redshift discs only, $z=5.6-4.0$. Indeed, we read $f\simeq 0.6$ at $T=3\times 10^{4}$K as in Fig. 3. There is a step near $10^{4}$K, resulting from the drop in the cooling curve near that temperature and from the fact that star formation is allowed below this temperature. The resultant plateau between $10^{4}$ and $10^{6}$K (the typical virial temperature), which is totally flat for the high-redshift galaxies, implies that when the threshold temperature is varied in this range the hot filling factor varies by only about $\pm 0.05$. The self-regulated value of $f$ is thus robust to the choice of the temperature threshold in the given range, with slightly smaller values of $f$ for higher $T$ thresholds.

Figure 5 shows the global 3D KS relation for all the discs of the VELA simulations, either using the whole disc, or using nine patches within each disc. The macroscopic gas density refers to cold gas, $T<3\times 10^{4}$K. The SFR is determined using the stars younger than $40-80\,{\rm Myr}$, in the minimum-bias way outlined in §A and in Tacchella et al. (2016b). In both panels we see a tight KS relation, with a slope $s\simeq 1.5$. The slope turns out to be rather insensitive to the way the SFR is computed. The colors, which refer to redshift, indicate that there is no significant variation of the KS relation with time, despite the systematic evolution from highly perturbed discs at high redshift to more relaxed discs toward $z\sim 1$. No explicit lower limit is applied to the gas density in the cells (beyond the upper limit to the temperature), motivated by the notion that the formation of molecular gas is not properly resolved on the grid-cell level. When a threshold of $n>1\,{\rm cm}^{-3}$ is applied, the slope of the KS relation becomes slightly flatter, $s\simeq 1.2-1.3$.

The fact that the KS slope turns out to be almost exactly $1.5$, as predicted in §2 based on the standard evolution of single SN bubbles, is not to be over-interpreted, as the exact evolution of SN bubbles in their early phases is not resolved in these simulations, and as additional radiative feedback is incorporated. Furthermore, one cannot (yet) reject the possibility that the macroscopic slope reflects the local SFR recipe as incorporated in these simulations, although the latter imposed a density threshold, which should have modified the power-law relation. The robustness of the macroscopic KS relation to the local SFR recipe should be explored via simulations where different local SFR recipes are incorporated (§4). There are preliminary indications from VELA cosmological simulations that the KS relation is robust to differences in the feedback recipes (Ceverino et al., 2014), which we will revisit in §4 using isolated-galaxy simulations. The meaningful new finding from the cosmological simulations is the self-regulation to a constant hot filling factor, Fig. 3, which, for a given density dependence of the filling factor, may be the main driver of the KS relation.

4 Isolated simulations: Role of SN Feedback

Robust to Local SFR

The validity of self-regulation to a constant hot filling factor, the role of SN feedback in it, and the insensitivity of the global KS relation to the local SFR recipe, are tested here via hydro-gravitational simulations of isolated galaxies.

4.1 The isolated-galaxy simulations

Our idealized simulations of isolated galaxies are carried out with the RAMSES (Teyssier, 2002) adaptive mesh refinement (AMR) code, including self-gravity, hydrodynamics, cooling and heating as well as star formation and stellar feedback based on subgrid recipes (Renaud et al., 2013). The two sets of simulations mimic star-forming disc galaxies at redshifts $z\sim 0$ and $z\sim 2$, using gas fractions of 15% and 50%, respectively. They are similar to those presented in Bournaud et al. (2014) but with several combinations of star formation and feedback recipes. Face-on and edge-on images of mass-weighted gas density along the line of sight in our fiducial simulation with 50% gas fraction are shown in Fig. 6 for a visual impression.

The simulation box size is $118\,{\rm kpc}$, the largest grid-cell size is $230\,{\rm pc}$ and the smallest cell after maximum refinement is $3.6\,{\rm pc}$. AMR cells are refined as soon as they contain more than 100 particles and/or a gas mass larger than $9\times 10^{3}M_{\odot}$ where the density is below $0.3\,{\rm cm}^{-3}$, or when they gas mass is larger than $2.8\times 10^{3}M_{\odot}$ with a density above $0.3\,{\rm cm}^{-3}$. Gas denser than $0.03\,{\rm cm}^{-3}$ is thus refined at $115\,{\rm pc}$ or better, and gas denser than $300\,{\rm cm}^{-3}$ is refined at the maximum, $3.6\,{\rm pc}$ resolution. The grid is also refined when the local Jeans length is smaller than 4 times the local cell size, and a density-dependent temperature floor keeps the Jeans length resolved by at least 4 cells at the highest refinement level (Teyssier, Chapon & Bournaud, 2010). Cooling is tabulated at solar metallicity and heating from a uniform UV background is included, as in Bournaud et al. (2014).

Each galaxy starts with a baryonic mass of $8\times 10^{10}M_{\odot}$, distributed in a stellar disc, a stellar bulge, and a gas disc. The gas mass is 15% or 50% of the baryonic mass. Each of the disc components has an exponential radial density profile with a scale length of $4\,{\rm kpc}$, truncated at $12\,{\rm kpc}$. The density profile vertical to the disc is exponential with a scale length of $600\,{\rm pc}$ for the stars, and $600\,{\rm pc}$ or $150\,{\rm pc}$ for the gas in the gas-rich and gas-poor cases, respectively. The initial gas temperature is $5\times 10^{4}$K. The bulge contains 20% of the stellar mass. It is spherical, obeying a Hernquist (1990) profile with a radial scale length of $600\,{\rm pc}$, truncated at $2\,{\rm kpc}$. Each galaxy is embedded in a spherical dark-matter halo with a Burkert (1995) density profile of a characteristic radius $20\,{\rm kpc}$. The mass of the dark matter within the effective disk radius is set to be half the baryon mass within that radius.

Local star formation is allowed in cells of hydrogen number density above a threshold $n=10\,{\rm cm}^{-3}$. Above this threshold, the SFR density is modeled as

[TABLE]

In the “standard” scheme, the recipe follows $\rho_{\rm sfr}=\epsilon_{\rm ff}\rho_{\rm g}/t_{\rm ff}$ (eq. 1), namely $s_{\rm loc}=1.5$. The normalization is determined by setting the efficiency to $\epsilon_{\rm ff}=0.025$. To test for sensitivity of the global slope to the local SFR recipe, the same simulations were run with a “steep” local slope of $s_{\rm loc}=2$ and with a “shallow” local slope of $s_{\rm loc}=1$. The normalization in each cases was tuned to obtain the same SFR as in the standard case during the first $20\,{\rm Myr}$.

Feedback is modeled through a combination of supernova and radiative feedback. The supernova feedback models the explosions of type-II supernovae, assumed to explode $10\,{\rm Myr}$ after the star birth. Seventy five percent of the initial energy of each supernova, assumed to be $10^{51}\,{\rm erg}$, is injected in the $\sim 6\,{\rm pc}$ vicinity of the SN in the form of thermal energy (50%) and kinetic energy (25%), while the rest is assumed to have radiated away before the shell reached $6\,{\rm pc}$. This is a typical energy budget for supernovae explosions once they have expanded to $\sim 6\,{\rm pc}$ (Martizzi, Faucher-Giguère & Quataert, 2015). The radiative feedback, termed HII feedback, models the photo-ionization of gas by O/B stars and the radiation pressure on the gas and dust. We assume a Strömgren sphere approximation around stellar particles younger than $10\,{\rm Myr}$ (see Renaud et al., 2013, for a more elaborate description). To test the role of feedback, and the contributions of SN feedback versus HII feedback, we ran the same simulations, with the standard local SFR recipe, but with four variants of feedback recipes; the standard SN+HII feedback, SN feedback alone, HII feedback alone, and no feedback. In all cases the simulation started with the fiducial feedback recipes, to allow relaxation from the initial conditions to a similar realistic configuration with spiral arms. Then, near $t=200\,{\rm Myr}$, if needed, the recipe was changed to the desired form.

The galaxy is allowed to relax for $\sim 200\,{\rm Myr}$ from the initial conditions to a realistic configuration, involving disc instability and spiral arms. The analysis is performed after this time for several hundred Megayears (typically until $580\,{\rm Myr}$), after which the approximation of an isolated galaxy gradually becomes less valid. The relevant quantities are measured in boxes that are spread throughout the disc, centered on the disc central plane. The box size is $1\,{\rm kpc}\times 1\,{\rm kpc}\times h$, where $h=0.5\,{\rm kpc}$ and $0.1\,{\rm kpc}$ for the high-$z$ and low-$z$ galaxies respectively. Alternatively we refer to the whole disc, within the disc exponential scale-length and full height $h$. In each box, we compute the SFR and supernova rate densities $\rho_{\rm sfr}$ and $S$, the cold gas number density $n$ for $T<5\times 10^{4}$K, and the volume filling factor for hot gas of $T>5\times 10^{4}$K. As in the cosmological simulations, the results are found to be insensitive to changing the temperature threshold by a factor of two.

4.2 Results: Role of SN Feedback and Robustness to Local SFR

Figure 7 shows $S$ versus $n$ (left) and the distribution of $f$ in the sub-volumes of each galaxy with the fiducial local SFR and feedback recipes, measured in four random snapshots in the time range $t=200-400\,{\rm Myr}$. A KS relation is demonstrated, with a tight distribution about $S\propto n^{1.5}$. In principle, this might be built in by the standard local SFR recipe assumed in the simulations. The non-trivial result is that the filling factor is distributed in a narrow range, $f\simeq 0.4\pm 0.07$, with only a small deviation between the low and high-$z$ galaxies. This is consistent with the idea that the system self-regulates itself to a constant filling factor of $\sim 0.5$, and points to the actual value at the self-regulated state for the isolated-galaxy simulations.

Figure 8 shows the time evolution of the hot filling factor in each variant of the two isolated galaxies (left and right columns). The solid curves refer to the standard SFR and feedback recipes. After an initial adjustment period of $\sim 200\,{\rm Myr}$, where the galaxy relaxes from the initial conditions, the filling factor becomes self-regulated to a rather constant value, similar in the two galaxies independent of the gas fraction. This self-regulation to a constant value is similar to what we saw in the cosmological simulations, though the asymptotic value of $f$ is somewhat lower for the isolated galaxies, reflecting the different SFR and feedback recipes and the other differences between the two simulation types performed with totally different codes (ART and RAMSES), initial conditions and environments.

With the feedback turned off in the bottom panels of Fig. 8, the hot filling factor becomes much lower. The non-zero value is likely due to shock heating within the supersonic turbulent ISM, and is possibly partly a remnant of the initial period with full feedback. This demonstrates that the hot filling factor predominantly represents hot bubbles generated by feedback, as assumed throughout this paper. SN-feedback alone gives rise to a hot filling factor similar to the case of full feedback, while HII-feedback alone leads to a lower filling factor, almost as low as with no feedback. This indicates that SN feedback dominates over HII feedback (§7).

The top panels of Fig. 8 demonstrate that the convergence of the hot filling factor to a constant value is robust to variations in the slope of the local SFR recipe, with small variations in the value of this constant. Together with the feedback dependence shown in the bottom panels, this indicates that the self regulation is robust and is driven by SN feedback, confirming the basic ansatz of our model. The typical hot filling factor at $t>200\,{\rm Myr}$ is listed for each run in Table 1.

Figure 9 and Fig. 10 present the global KS relations as produced in the isolated-galaxy simulations with gas fractions 0.15 and 0.50 respectively. The SFR density and cold-gas density are measured in macroscopic volumes within the disc at different snapshots after the initial $225\,{\rm Myr}$. Table 1 lists the slope of the global KS relation in the different cases, as determined by a linear fit for $n>50\,{\rm cm}^{-3}$, and the rms scatter about the linear relation. With the standard local SFR recipe and feedback (top-left, black), the global KS relation has a slope $s\simeq 1.43$, with a small scatter of $\pm 0.06$dex. When the local slope is steeper ($s_{\rm loc}=2$) or flatter ($s_{\rm loc}=1$), the top panels show that the global slope becomes only slightly steeper or flatter, by $\simeq 0.1$, with a scatter $\pm 0.13$dex. In general, both for the simulations with low and high gas fraction, the global slope remains roughly the same and with a small scatter, insensitive to the local slope.

The bottom panels of Fig. 9 and Fig. 10 show that, for the standard local SFR recipe, the KS relation is almost the same with the full feedback and with SN feedback alone. When only HII feedback is activated, the global slope steepens to $s\simeq 1.88$, and the scatter grows to $\pm 0.15$dex. With no feedback, the global slope steepens further to $s\simeq 2$, and the scatter grows further to $\pm 0.21$dex. This is consistent with our basic assertion that SN feedback dominates over the HII feedback, and is responsible for the KS relation (see §7).

The weak dependence of the large-scale KS relation on the small-scale SFR recipe is consistent with earlier tests using SPH hydro simulations of isolated galaxies (Hopkins, Quataert & Murray, 2011; Hopkins, Narayanan & Murray, 2013). Considering galaxies of a variety of masses at low and high redshifts, they experimented with a variety of local star-formation recipes. Hopkins, Quataert & Murray (2011) explored density criteria with a SFR efficiency in the range $\epsilon_{\rm ff}=0.0035-0.06$, a power-law dependence on gas density in the range $\rho_{\rm sfr}\propto n^{1.0-2.0}$, and a threshold gas density for star formation in the range $n_{\rm min}=10-2500\,{\rm cm}^{-3}$. Hopkins, Narayanan & Murray (2013) experimented with different local physical criteria for star formation, including self-gravity, Jeans instability, density, temperature, molecular-gas content and cooling rate. They found that once feedback is incorporated, the galaxy is self-regulated to a global KS relation with only little sensitivity to the local SFR recipe. This is in the sense that for a galaxy of given global galaxy properties, such as gas density, the global SFR approaches roughly the same value independent of the local SFR recipe. Hopkins (2014) then showed that in their FIRE cosmological simulations, using their fiducial SFR recipe, the simulation converges to a KS relation similar to the observed relation. These experiments, like ours, indicate that the SFR is robustly regulated by feedback. In the current paper, we address the origin of this through the robust self-regulation of the bubble hot volume filling factor into a constant value. Our simulations demonstrate the convergence to a constant filling factor both in cosmological and isolated settings, and our isolated-disc simulations confirm the robustness to different local SFR recipes, and establish the dominant role of SN feedback in the self-regulation to a constant filling factor and the generation of the global KS relation.

5 Clustered Supernovae - Analytic

In the previous two sections we established the validity of the concept of self-regulation into a constant hot filling factor in galaxy simulations that incorporate star-formation and SN feedback in a realistic ISM, with the resultant KS relation independent of the local star-formation recipe and apparently driven by SN feedback. We now return to simplified analytic modeling, making first steps in generalizing the analytic modeling of isolated SNe (§2). Here we address idealized cases of clustered SNe in star-forming clouds, where a sequence of SNe explode at the same location. In the next section we test these analytic estimates with spherical hydrodynamical simulations. Then, in the following section, we consider the Strömgren bubbles photo-ionized by the pre-SN O/B stars and their interplay with the SN bubbles. These simplified analytic attempts will help verifying the validity of the basic model for the origin of the KS relation in somewhat more realistic circumstances, and physically interpreting the behavior in the full simulations.

5.1 Four characteristic times

As star formation tends to occur in clumps (molecular clouds, star clusters), many SNe occur practically in the same place. This is likely to affect the bubble filling factor and it may change its density dependence and therefore the resultant KS relation. We attempt here to estimate the possible effects of clustering in a simplistic analytic way, to be followed by simplified simulations. We model the sequence of SNe in a cluster with two independent parameters. We consider first a situation where $\nu$ SNe occur in each point-like cluster, with a constant SN rate during a burst duration $t_{\rm b}$, associated with the lifetime of the star-forming cluster. The typical time available for a SN before the successive SN explodes is thus $t_{\rm s}=t_{\rm b}/\nu$, so the two parameters could be $t_{\rm s}$ and $t_{\rm b}$. For a given cluster of SNe we consider the bubble about it. In analogy to the case of an individual SN bubble, the cumulative bubble will have a phase analogous to the adiabatic phase until a cooling time $t_{\rm c}$, followed by phases analogous to the snow-plow phase in which the outer shell has collapsed to a thin massive shell pushed by a wind or pressure, whose speed eventually fades away to the ISM sound speed at a fading time $t_{\rm f}$.

The evolution of the cumulative bubble, its fading time and radius, the resulting bubble filling factor, and the final power index $s$ of the density dependence in the KS relation, $\rho_{\rm sfr}\propto n^{s}$, depend on the interplay between the timescales characterizing the SN cluster, $t_{\rm s}$ and $t_{\rm b}$, and those of the cumulative bubble, $t_{\rm c}$ and $t_{\rm f}$. Note that by definition $t_{\rm s}<t_{\rm b}$ and $t_{\rm c}<t_{\rm f}$.

5.2 Different zones in parameter space

We divide the parameter space into the following different zones, where the analysis of the bubble consists of different phases such that the final density dependence may be different (see a schematic cartoon in Fig. 11):

A. A short burst, $t_{\rm b}<t_{\rm c}$

A1. instantaneous, $t_{\rm b}\ll t_{\rm c}$

A2. non-instantaneous, t_{\rm b}\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}t_{\rm c}

B. **Long and continuous, $t_{\rm s}<t_{\rm c}<t_{\rm b}$

**B1. moderately long, $t_{\rm b}<t_{\rm f}$

B2. very long, $t_{\rm f}<t_{\rm b}$

C. **Long and semi-continuous, $t_{\rm c}<t_{\rm s}<t_{\rm f}$

**C1. moderately long, $t_{\rm b}<t_{\rm f}$

C2. very long, $t_{\rm f}<t_{\rm b}$

D. Very long and separable, $t_{\rm f}<t_{\rm s}$.

We summarize here the expectations in each zone, and elaborate in the following subsections.

For a short burst, zone A, the solution can be deduced from the solution for an individual SN, namely a KS relation with a power index $s=1.48$. This is obvious for a very short burst, case A1, where the SN cluster is analogous to a single hyper-nova with the energy multiplied by $\nu$. We show below that the same power law with $s=1.48$ is expected also in the non-instantaneous zone A2.

A similar power index of $s=1.48$ is expected also in the opposite extreme of a very long time interval between SNe, zone D. If each individual bubble fades away well before the explosion of the subsequent SN, such that the faded bubble had time to recover the original unperturbed ISM environment, the cluster is expected to behave like a sequence of separable individual SNe. The solution is the single-SN solution, with a power index $s=1.48$ in the KS relation. This may still be a sensible crude approximation when $t_{\rm s}\sim t_{\rm f}$, near the border of zones D and C.

In the sub-cases of zone B, we predict analytically slopes in the range $s=1.1\!-\!2.0$, and typically find values closer to $s=1.5$ in the simulations. Here, as well as in zone A2, as long as $t<t_{\rm c}$, and actually as long as $t<t_{\rm b}$, the energy source by the frequent SNe can be treated as a continuous energy flux, causing a wind-driven expansion phase at a rate $R_{\rm s}\propto t^{3/5}$. In zone B2, where the fading occurs while the burst is still on, $t_{\rm f}<t_{\rm b}$, our analytic estimate leads to $s=2.0$. In zone B1, where the fading occurs well after the burst is over, $t_{\rm b}<t_{\rm f}$, after the bubble went into a passive snow-plow-like phase of $R_{\rm s}\propto t^{2/7}$ at $t>t_{\rm b}$, our analytic estimate becomes $s=1.1$. In the intermediate situation between B1 and B2, where $t_{\rm f}\sim t_{\rm b}$, we expect a value of $s$ between $1.1$ and $2.0$, growing as a function of $t_{\rm b}/t_{\rm f}$.

In zone C, the continuous limit may still be valid with certain modifications, and the trend of the slope $s$ with $t_{\rm b}/t_{\rm f}$, from C2 to C1, is likely to be similar to the trend in cases B.

We describe below these analytic predictions, and follow in §6 with spherical simulations that confirm the convergence to $s\sim 1.5$ in most of the different cases.

5.3 A Continuous Energy Source

Once $t_{\rm s}<t_{\rm c}$ (zones A and B), and possibly even as long as $t_{\rm s}<t_{\rm f}$, (zone C), while the burst is on, $t<t_{\rm b}$, we follow Weaver et al. (1977) in treating the energy source from the successive SNe as continuous, with a constant power (“luminosity”),

[TABLE]

Here $\gamma$ is the efficiency by which each SN delivers energy $E$ in the clustered environment, assumed to be of order unity. Hereafter, we let $e_{51}$ actually represent $\gamma\,e_{51}$. We can write

[TABLE]

where $L=10^{38}\,{\rm erg}\,{\rm s}^{-1}\ell_{38}$ and $t_{\rm s}=1\,{\rm Myr}\,{t_{\rm s}}_{6}$.

The single-SN Sedov-Taylor expansion, $R_{\rm s}\propto t^{2/5}$, is now replaced by a slightly faster expansion, $R_{\rm s}\propto t^{3/5}$, as the constant $E$ in eq. (6) is replaced by $E(t)=L\,t$ where $L$ is a constant. The shell properties prior to the cooling time become222The reduced factor of 0.88 in eq. (30) compared to the single SN case is because in the case of a shock driven by a constant wind only 55% of the energy is made available to the swept-up ISM gas while the rest is stored in the shocked wind behind it. The thermal pressure of this shocked wind is pushing the swept-up mass, and this is somewhat less efficient than in the Sedov-like expansion.

[TABLE]

[TABLE]

[TABLE]

where $t=1\,{\rm Myr}\,t_{6}$.

Proceeding in analogy to the single-SN case, eq. (10), while considering the total energy injected so far, $E(t)=Lt$, the cooling time and the shell radius and velocity at cooling become

[TABLE]

[TABLE]

[TABLE]

where the spatial integral in eq. (10) has been evaluated numerically for the wind-like similarity solution. At this time the shell exits the simple adiabatic continuous regime, where it expands $\propto t^{3/5}$ according to eq. (30), and it may enter another regime, depending on the zone in parameter space.

5.4 Zone A: Short Burst $t_{\rm b}<t_{\rm c}$

For an instantaneous burst, $t_{\rm b}\ll t_{\rm c}$, the cluster is analogous to a single SN with an energy $\nu E$. For a given overall SN rate $S$, the rate of clustered explosions is $S/\nu$. According to the dependence of the filling factor on energy in eq. (21), the filling factor in the clustered case is

[TABLE]

Thus, the KS relation obtained by requiring $f_{0}=const.$ is similar to that of single SNe, eq. (24), with the same density dependence $\rho_{\rm sfr}\propto n^{1.48}$, and with the amplitude, or $\epsilon_{\rm ff}$, scaled as $\nu^{-0.26}e_{51}^{-1.26}$. For example, with GMCs of $10^{4}M_{\odot}$, and $\mu_{\rm sn}=100\,M_{\odot}$ stellar mass per supernovae, one has $\nu=100$, so the supernova rate at a constant $f_{0}$, and similarly $\epsilon_{\rm ff}$, is smaller by a factor $\nu^{0.26}\simeq 3.3$ compared to the case of unclustered SNe (assuming $\gamma=1$).

In general in zone A, $t_{\rm b}<t_{\rm c}$, as long as the burst is on, $t<t_{\rm b}$, the shock is driven by a continuous-energy wind, $R_{\rm s}\propto t^{3/5}$, following eq. (30).

Once the energy input ceases but before the shell losses a significant fraction of its energy, $t_{\rm b}<t<t_{\rm c}$, the shock enters a Sedov-Taylor-like expansion,

[TABLE]

[TABLE]

where eq. (29) has been used to express $L$ and $t_{\rm b}$ by $\nu$. Note that while the numerical values are different, the dependence of these expressions on the variables are the same as the expressions for the single-SN Sedov-Taylor solution, eq. (6) and eq. (7), once the single-SN energy $E$ is replaced by the total energy $\nu\gamma E$.

In computing the cooling time, the energy loss to radiative cooling is the sum of the integrals of $\dot{E}[R_{\rm s}(t)]$, eq. (10), over time in the successive intervals $0-t_{\rm b}$ and $t_{\rm b}-t$, where the wind phase and the Sedov-Taylor-like phase are assumed to be valid respectively. In the spatial integrals of eq. (10) in each phase, the density and temperature profiles are obtained from the corresponding self-similar solutions, and the upper bound for the integration, $R_{\rm s}(t)$, is from eq. (37) and eq. (6) respectively. After some algebra, The cooling time, where the energy loss is one third of $E(t)$, turns out to be

[TABLE]

[TABLE]

In the limit $t_{\rm b}\ll t_{\rm c}$, zone A1, the second term is negligible, and one recovers the single-SN case with the energy multiplied by $\nu\gamma$, as expected. In the other limit of zone A, $t_{\rm b}\sim t_{\rm c}$, an estimate is obtained by substituting $t_{\rm b}=t_{\rm c}$,

[TABLE]

After the cooling time, at $t>t_{\rm c}$, the shell is in a passive snow-plow phase,

[TABLE]

[TABLE]

The shell fades away when the shock velocity is reduced to the sound speed of the medium,

[TABLE]

[TABLE]

With a rate of $S/\nu$ for the clusters of SNe, the “hot” volume filling factor becomes

[TABLE]

Recall that some of the $n$ dependence is in $K$, with the same sign as the explicit $n$ dependence in eq. (46). The resultant power of the $n$ dependence in zone A2 is thus always $s>1.12$. When the second term in eq. (40) is small, the power is $s\rightarrow 1.47$. The cubic power of $t_{\rm b}$ in $K$ makes this a good approximation almost all the way to $t_{\rm b}\sim t_{\rm c}$. We conclude that throughout zone A the KS relation is reproduced with $s\simeq 1.5$, and with the normalization scaling as $\nu^{-0.26}$.

5.5 Zone B: Long, Continuous Burst $t_{\rm s}<t_{\rm c}<t_{\rm b}$

As long as $t<t_{\rm c}$, the shock is driven by a continuous-energy wind, $R_{\rm s}\propto t^{3/5}$, following eq. (30). The cooling time, radius and velocity are given by eq. (33), eq. (34) and eq. (35).

After the cooling time, during $t_{\rm c}<t<t_{\rm b}$, the shocked ISM has cooled, lost its pressure, and collapsed to a thin shell, which is now pushed outward by the pressure of the shocked wind region. Radiation losses in this region can be ignored as the typical gas velocity is $\sim 2000\,{\rm km}\,{\rm s}^{-1}$ and the temperature is $5\times 10^{7}$K. The total energy of the shocked wind region is (Weaver et al., 1977, eq. 14)

[TABLE]

Since this energy is related to the pressure via

[TABLE]

the pressure is

[TABLE]

The energy change in the shocked wind region includes the work done by the expanding shell,

[TABLE]

The shell equation of motion is then

[TABLE]

eqs. (49), (50) and (51) can be solved to obtain the shock radius and velocity

[TABLE]

[TABLE]

This is similar to eq. (30) and eq. (31), which describe the evolution of the shell in the earlier phase driven by a continuous wind before cooling, at $t<t_{\rm c}$. The pre-factors are smaller here due to the collapse of the cooling shell, during which the shock temporarily slows down.

In order to address the fading one should consider two different cases where the fading occurs either before or after the end of the burst at $t_{\rm b}$.

5.5.1 Zone B2: Very long burst $t_{\rm f}<t_{\rm b}$

If the burst is very long, the shock velocity could reach the ISM speed of sound while the burst is still on, during the active snow-plow phase. In this case, based on eq. (52) and eq. (53), the fading would occur at

[TABLE]

[TABLE]

The “hot” filling factor is therefore

[TABLE]

Thus, in zone B2, under the approximations made, the predicted slope for the KS relation is $s=2$. The normalization of th KS relation now scales as $\nu^{-1}$, but also as $t_{\rm b}^{2}$, and with a very strong dependence on the sound speed $c^{7}$.

5.5.2 Zone B1: Moderately long burst $t_{\rm b}<t_{\rm f}$

If the burst is not so long, such that the fading occurs only after the burst is over, $t_{\rm b}<t_{\rm f}$, the pre-faded shell eventually enters a regime, $t>t_{\rm b}$, where its evolution is no longer governed by eq. (52). In this regime, the dynamics of the shell is similar to the passive snow-plow phase of a single SN,333The switching off of the driving of the shell and the transition to a snow-plow similar to a single SN can be considered to be instantaneous with respect to $t_{\rm b}$ since the information about the turning off of the energy source is transferred from the center to the shell on a short timescale of $\sim 0.1\,{\rm Myr}$, traveling a distance $R_{\rm s}\sim 100\,{\rm pc}$ at a speed $c_{\rm s}\sim 1000\,{\rm km}\,{\rm s}^{-1}$.

[TABLE]

[TABLE]

The fading, $V_{\rm s}=c_{\rm s}$, thus occurs at

[TABLE]

[TABLE]

where $\ell_{38}$ was inserted using eq. (29). The “hot” volume filling factor is therefore

[TABLE]

Hence, in zone B1, under the approximations made, the predicted slope for the KS relation is $s=1.12$. The normalization scales weakly with $\nu^{-0.12}$ and more significantly with $t_{\rm b}^{-0.64}$ and $c^{2.6}$.

We have not worked out analytic predictions for zone C, where $t_{\rm c}<t_{\rm s}<t_{\rm f}$. The spherical simulations described in §6 indicate that the behavior in zone C is similar to the behavior in zone B (B1 or B2 respectively). One limitation here is that for a given cluster of SNe the bubble may be in different zones of parameter space for different $n$ values, as the evolution depends on $t_{\rm s}/t_{\rm f}$, and $t_{\rm f}$ varies with $n$.

5.6 More Realistic Clusters of SNe

Before proceeding to simulations that may refine the analytic estimates of this section, it is worth evaluating the possible assignment of actual observed star-forming clusters to zones in parameter space according to the above scheme.

5.6.1 A uniform time sequence

First consider a uniform sequence of SNe, as estimated so far, which may serve as a crude approximation for real clouds in some cases. For example, we estimated for the Milky Way in eq. (12) and eq. (18) that for a single SN $t_{\rm c}\sim 0.05\,{\rm Myr}$ and $t_{\rm f}\sim 4\,{\rm Myr}$. The giant molecular cloud (GMC) lifetimes are typically tens of Myr, after which they are disrupted by feedback, so we may adopt $t_{\rm b}=50\,{\rm Myr}\,{t_{\rm b}}_{50}$. This puts the GMCs in the “long-burst” zones, B to D, where $t_{\rm c}<t_{\rm b}$.

In order do further distinguish between zones B to D, we should estimate the time between SNe $t_{\rm s}$. Considering a cluster stellar mass $M_{\rm c}=10^{4}M_{\odot}\,M_{4}$, and a stellar mass per SN $\mu_{\rm sn}=100\,\mu_{2}$, the number of SNe is $\nu=100\,M_{4}\mu_{2}^{-1}$, which gives

[TABLE]

If the cloud lifetime $t_{\rm b}$ is much shorter than $50\,{\rm Myr}$, or for massive clouds of $\sim\!10^{5}M_{\odot}$, one may have t_{\rm s}\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}t_{\rm c}, namely the cloud is in zone B. Otherwise, for clouds that live longer or are less massive, one has $t_{\rm c}<t_{\rm s}$, namely the cloud is in zone C. In extreme cases, of long-lived, low-mass GMCs, or if $\mu_{\rm sn}$ is somehow particularly large within GMCs, the time per SN, $t_{\rm s}$, could be comparable to $t_{\rm f}$ such that we are barely in zone D. We thus expect the low-mass, intermediate and massive GMCs, under the uniform-sequence approximation, to be in zone D, C and B respectively. If $t_{\rm f}<t_{\rm b}$ they would likely be assigned to zone B2 or C2.

In high-redshift giant clumps the lifetime is expected to be longer, possibly $t_{\rm b}\sim 300\,{\rm Myr}$, the characteristic migration time within the violently unstable discs into the central bulge (Dekel, Sari & Ceverino, 2009). On the other hand the clump masses are expected to be much larger, comparable to the Toomre mass, e.g., $M_{\rm c}\sim 10^{8}M_{\odot}$ (Dekel, Sari & Ceverino, 2009; Mandelker et al., 2017). This gives $t_{\rm s}\sim 3\times 10^{-4}\,{\rm Myr}$, namely $t_{\rm s}\ll t_{\rm c}$, so the clump is in zone B. With $t_{\rm f}<t_{\rm b}$, the high-$z$ giant clumps are likely to be assigned to zone B2.

5.6.2 A Cloud of Clusters

To be more realistic, we wish to evaluate the evolution of the cumulative super-bubble about a star-forming cloud of gas (g) that consists of smaller star clusters (c), each generating a short burst. We assume the number of clusters in the cloud to be $\nu_{\rm g}=10\,\nu_{{\rm g},10}$. For a cloud of mass $M_{\rm g}=10^{5}M_{\odot}\,M_{5}$, assuming one SN per $100M_{\odot}\,\mu_{2}$ (typically $\mu_{2}\sim 1.5$), we expect a total of $\nu=10^{3}M_{5}\mu_{2}^{-1}$ SNe in the cloud, namely $\nu_{\rm c}=100M_{5}\mu_{2}^{-1}\nu_{{\rm g},10}^{-1}$ SNe per cluster.

Within a cluster, we write the duration of each burst of SNe as $t_{\rm b,c}=1\,{\rm Myr}\,t_{\rm b,c,6}$. Dividing by $\nu_{\rm c}$, the average time between successive SNe is

[TABLE]

Recall that the individual SN cooling time, from eq. (12), is

[TABLE]

and the individual fading time, from eq. (18), is

[TABLE]

We learn that, within each cluster,

[TABLE]

This indicates that the cluster is not in zone A.

To test whether the cluster could be in zone D, we examine the ratio

[TABLE]

This implies that the cluster is not in zone D unless the burst is very long, the cloud is of much lower mass than $10^{5}M_{\odot}$, and the hydrogen density is very high.

To distinguish between zones B and C, we examine the ratio

[TABLE]

This implies that for massive clouds and relatively short bursts the individual clusters would be in zone B. If the SN burst is much longer than $\sim 1\,{\rm Myr}$, or the cloud is significantly less massive than $\sim 10^{5}M_{\odot}$, each cluster would be in zone C.

If the cluster is in zone B, in order to distinguish between zones B1 and B2, we insert $\nu_{\rm c}$ in the expression for $t_{\rm f}$ in zone B1, eq. (59), and obtain

[TABLE]

This is zone B1 if $t_{\rm b,c,6}<10$ and all other factors are unity. More so if $c_{1}<10\,{\rm km}\,{\rm s}^{-1}$ and if $M_{5}>1$. However, it may be zone B2 if $t_{\rm b,c,6}>10$, $M_{5}\ll 1$, $\gamma\ll 1$, $\nu_{{\rm g},10}\gg 1$ or $n_{0}\gg 1$.

In the cloud of clusters, we assume that the duration of the burst of clusters is $t_{\rm b,g}=50\,{\rm Myr}\,t_{\rm b,g,1.7}$. This implies a long average duration between successive clusters of

[TABLE]

To test whether the cloud is in zone D, we compare $t_{\rm s}$ with the cumulative fading time as evaluated for each cluster in zone B1, from eq. (59) with $\nu_{\rm c}$, and obtain

[TABLE]

This is of order unity, and can be larger (zone D) if $M_{5}$ is low and if $n$ is high.

6 Clustered Supernovae - Spherical Simulations

6.1 Method of Spherical Simulations

We performed 1D spheri-symmetric hydrodynamical simulations about a fixed center to test the idealized analytic estimates and to extend the results to cases where analytic modeling is more difficult. The simulations utilize the finite volume Eulerian hydrodynamical code pluto-v4.0 (Mignone et al., 2007), which solves the standard conservation equations for mass, momentum and energy, with gravity turned off. The energy is kinetic and thermal, with $\gamma=5/3$, and it changes by input from the SN and radiative losses. The background density is kept uniform, at $T=10^{4}$K, with a sound speed of $c_{\rm s}=[\gamma k_{\rm b}T/(\mu m_{\rm p})]^{1/2}\simeq 15.1\,{\rm km}\,{\rm s}^{-1}$.444This is assuming $\gamma=5/3$ and $\mu=0.6$, which is valid for T\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}2\times 10^{4}K, but a constant value of $\mu$ is adopted for simplicity. The simulations are stopped once the shock speed is reduced to $c_{\rm s}$. The metallicity is kept at the Solar value. The grid resolution is uniform throughout the box and constant in time. The default resolution is $\Delta r=0.0244$ and $0.244\,{\rm pc}$ for single SNe and clusters respectively, and it is increased by factors up to 20 to test for convergence (see Table 2).

A SN is assumed to inject in the ISM gas mass of $M_{\rm sn}=5M_{\odot}$ and thermal energy of $E_{\rm sn}=10^{51}\,{\rm erg}$, at a rate that is uniform within a spherical volume of radius $r_{\rm inj}\sim 1\,{\rm pc}$ about the center and during a short time interval $\Delta t_{\rm inj}$, selected to be the largest between the hydro timestep and $1\,{\rm kyr}$. In the clustered SNe runs, the first SN energy is injected in a region of $r_{\rm inj}=0.5\,{\rm pc}$, but the subsequent SNe are put in a larger region of $r_{\rm inj}=2\,{\rm pc}$ to avoid numerical instabilities when injecting energy in a very low density medium. The results are not sensitive to the exact value of $r_{\rm inj}$ and $t_{\rm inj}$ as long as they are much smaller than the cooling radius and time.555The small radius and time are chosen to avoid numerical effects that could arise if the SN energy was injected in a larger volume, such as a failure to generate a strong shock and/or having the injected gas at a temperature where it cools faster than it can expand. In order to avoid numerical instabilities due to sharp gradients of the thermodynamic quantities at $r_{\rm inj}$, the mass and energy injected as a function or $r$ are linearly smoothed in the shell $(0.9-1.1)r_{\rm inj}$. Radiative cooling is incorporated as in eq. (10) and eq. (9), assuming a solar metallicity.

6.2 Results of Spherical Simulations

6.2.1 Single SN

Figure 12 shows simulation results for a single SN, to test the success of the simulations in reproducing the textbook analytic predictions. The fiducial resolution is $\Delta r=0.024\,{\rm pc}$ and the energy injection is within $r_{\rm inj}=0.5\,{\rm pc}$. The left panel, which shows the evolution of the shock radius $R_{\rm s}(t)$ for $n=1\,{\rm cm}^{-3}$, demonstrates that the simulation recovers quite successfully the expected evolution in the adiabatic and snow-plow phases. The right panel shows the SN rate density $S$, assuming a constant bubble filling factor $f_{0}=0.6$, as a function of gas density $n$ for simulations with six different values of $n$. The filling factor is computed by $f_{0}=(4\pi/3)R_{\rm f}^{3}t_{\rm f}S$, where $t_{\rm f}$ and $R_{\rm f}$ are determined when the shock velocity have reduced to $15.1\,{\rm km}\,{\rm s}^{-1}$, the sound speed of the ISM. Then, assuming $f_{0}=const.$, $S(n)$ is determined from $R_{\rm f}(n)$ and $t_{\rm f}(n)$. In addition to the shock radius that marks the outer radius of the shell, we also determine the inner radius of the shell, $R_{\rm cd}$, encompassing the low-density hot bubble, defined where the shell density falls below half the ISM density $0.5\,n$. The results for the fiducial resolution are compared to the results from simulations with twice and twenty times better resolution, indicating convergence. We see that the simulations roughly reproduce the predicted slope of $S\propto n^{1.48}$, both for $R_{\rm s}$ and $R_{\rm cd}$, as well as the predicted normalization. This indicates that the 1D simulations could be useful for studying the SN bubbles in the clustered-SN cases.

6.2.2 Cases A

Figure 13 show simulation results in the $S-n$ plane, assuming $f_{0}=0.6$, for short clustered bursts in zones A1 and A2 of parameter space, where $t_{\rm b}\ll t_{\rm c}$ and t_{\rm b}\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}t_{\rm c} respectively. We simulate $\nu=10$ clustered SNe in each case. In the A1 case the burst duration is $t_{\rm b}=10^{-3}\,{\rm Myr}$, and in A2 it is chosen for each $n$ to be $0.5t_{\rm c}$, given that $t_{\rm c}$ depends on $n$, eq. (12). We could have also changed $\nu$ as we change $n$ in order to keep a constant luminosity, but the dependence of $f_{0}\propto\nu^{0.6}$ is rather weak, eq. (46). The fiducial resolution is $\Delta r=0.244\,{\rm pc}$ and the energy injection is within $r_{\rm inj}=2\,{\rm pc}$. The simulation results converge as the resolution is increased, and they reproduce the analytic predictions of §5.4, both in terms of slope, $S\propto n^{1.5}$, and amplitude.

6.2.3 Case B1

The simulation results for a clustered SNe case of type B1 are shown in Fig. 14, with $\nu=100$ SNe within $t_{\rm b}=1\,{\rm Myr}$. The fiducial resolution is again $\Delta r=0.244\,{\rm pc}$ and $r_{\rm inj}=2\,{\rm pc}$. Recall that the analytic model of §5.5.2 predicts an adiabatic wind phase of $R_{\rm s}\propto t^{3/5}$ until $t_{\rm c}$, followed by an active snow-plow phase of $R_{\rm s}\propto t^{3/5}$ driven by a continuous wind as long as the energy source is on, and a transition near $t_{\rm b}$ to a passive snow-plow with $R_{\rm s}\propto t^{2/7}$ until fading. At $t\sim t_{\rm s}$, the simulated growth of $R_{\rm s}$ is closer to the single-SN solution $t^{2/5}$ because of the discreteness of the few first SN explosions. The simulated growth then steepens into the expected continuous wind phase, but it shows a slightly slower growth rate during the active snow-plow after $t_{\rm c}$666This may partly arise from the assumption that the background pressure vanishes and partly be a numerical effect (Weaver et al., 1977) that is hard to overcome both in Eulerian and Lagrangian codes. and a somewhat steeper growth rate during the passive snow-plow, $R_{\rm s}\propto t^{0.42}$ compared to the expected $\propto t^{2/7}$, eventually reaching fading at roughly the same radius as predicted. A similar discrepancy between the simulations and the analytic model in the passive snow-plow phase is also seen in the Lagrangian simulations of Gentry et al. (2017), so it is likely due to an over-simplification in the model, which is yet to be understood. The right panel of Fig. 14 shows that for a constant filing factor, $f_{0}=0.6$, the simulations yield $S\propto n^{1.5}$, similar to the single-SN case and steeper than the $s=1.12$ predicted by the simplified analytic model of §5.5.2.

6.2.4 Case C2

It is difficult to predict analytically the evolution of $R_{\rm s}$ in zone C where $t_{\rm c}<t_{\rm s}<t_{\rm f}$ as neither a continuous wind solution can be applied (valid for $t_{\rm s}<t_{\rm c}$) nor the SNe can be considered as separable (valid for $t_{\rm f}<t_{\rm s}$). The evolution of $R_{\rm s}$ in a simulation in zone C is presented in Fig. 15, for $n=10\,{\rm cm}^{-3}$, with $t_{\rm c}\sim 0.1\,{\rm Myr}$, and $\nu=100$ SNe. The burst duration is chosen to be $t_{\rm b}=\xi\nu t_{\rm f}$ where $\xi=t_{\rm s}/t_{\rm f}$ varies from $1/3$ to $3/4$ in order to explore the behavior at t_{\rm s}\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}t_{\rm f}. The simulations show, as expected, that $R_{\rm s}$ evolves like a single SN, through adiabatic and passive snow-plow phases, till the second SN is injected. We learn that after a few more SNe have been injected, the expansion rate gradually steepens to the $R_{\rm s}\propto t^{3/5}$ expected for a continuous wind, as in zone B. Since we expect that for $t_{\rm f}<t_{\rm s}$ (where $t_{\rm f}$ is the single-SN fading time) the cluster would be in zone D, where the SNe can be treated as separate, we examine the behavior in zone C as a function of proximity of $t_{\rm s}$ to $t_{\rm f}$. We see in the figure that once $t_{\rm s}>0.7t_{\rm f}$ the shock stops expanding and it fades away (with a density jump of only 20%) after seven SNe (t\lower 2.58333pt\hbox{;\buildrel>\over{\sim};}7t_{\rm s}) despite the fact that SNe continue to explode at the centre. This implies that for $t_{\rm s}>0.7t_{\rm f}$ the cluster is already in the separable zone D, where the toy model predicts $s=1.5$.

It is difficult to study the $S-n$ relation through simulations in zone C, as we did in zone B, because the behavior is expected to depend on $t_{\rm s}/t_{\rm f}$ while $t_{\rm f}$ depends on $n$. Therefore, for a fixed $t_{\rm s}$ (given $t_{\rm b}$ and $\nu$, namely the same $\ell_{38}$), the bubble may be assigned to a different parameter zone for different $n$ values.

We note that a similar conclusion regarding the validity of the wind solution has been obtained by Gentry et al. (2017), using a Lagrangian code, who found that for clusters with $\nu\sim 10$ the wind solution in eq. (52) is not valid as t_{\rm s}/t_{\rm f}\lower 2.58333pt\hbox{;\buildrel>\over{\sim};}1, while an approximate wind solution prevails when \nu\lower 2.58333pt\hbox{;\buildrel>\over{\sim};}100.

6.2.5 A cloud of clusters

When simulating a cloud of clusters, we consider a total of 660 SNe (corresponding to a $10^{5}M_{\odot}$ super star cluster) over a time period of $30\,{\rm Myr}$. The SNe are divided into $\nu_{\rm g}=10$ clusters, each containing $\nu_{\rm c}=66$ SNe. The mass and energy of SNe inside a cluster are injected at time intervals of $t_{\rm s}=1/66\,{\rm Myr}$, while subsequent clusters are switched on at times separated by $3\,{\rm Myr}$, each lasting for $1\,{\rm Myr}$.

Figure 16 shows the $S-n$ relation, for a fixed $f_{0}=0.6$, from simulations of a cloud of clusters. We find the result to be, again, a good match to $S\propto n^{1.5}$. Eq. (69) tells us that for $n>3\,{\rm cm}^{-3}$ the fading time for the individual clusters is $t_{\rm f}<3\,{\rm Myr}$ (considering $t_{{\rm b,c},6}=1$ and $c=15.1$ km s*-1*), which is smaller than the separation between two consecutive clusters. This allows the clusters to contribute to the volume filling factor independently following a case B1 solution, which according to the simulation in Fig. 14 is $s\simeq 1.5$

6.3 Convergence of the simulations

Figure 12, referring to a single SN, seems to show in the left panel that for $n_{0}=1$, with the default resolution of $\Delta r=0.244\,{\rm pc}$, the simulated evolution of the shock radius recovers the theoretical prediction quite accurately. However, the right panel demonstrates that this is not the case at higher densities, where the simulated $S$ (namely the measured filling factor) at the default resolution (large filled light-blue symbols, sometimes hidden behind the blue symbols that are for a resolution twice as good) overestimates the theoretical prediction by $\sim 0.2$dex at $n_{0}=2$. This is a resolution effect, as at high ISM densities the shell density is also high, requiring higher resolution for a proper treatment of the cooling at the shock front and at the interface between the bubble and the shell. Increasing the resolution by a factor of 20 (blue-white symbols) makes the simulated $S$ approximate the theoretical solution to better than $0.1$dex even at $n_{0}=300$. Figure 13, left panel, where we also have an accurate theoretical model, shows a similar improvement of accuracy with improved resolution.

We thus conclude that our simulations for a single SN converge to the theoretical solution, and adopt the same resolution in the simulations of clustered SNe. It is interesting to note that the slope $s$ in the $S-n$ (namely in the KS) relation, as derived from the clustered-SN simulations, shown in Fig. 13 (right), Fig. 14 (right), and Fig. 16, at $n_{0}\geq 1$, is rather insensitive to the resolution, and can thus be derived from the simulations with the default resolution as well as from the simulation with higher resolution. We verified that our Eulerian simulations and the Lagrangian simulations of Gentry et al. (2017), when performed on a similar zone of parameters, roughly agree on the range of validity of the continuous wind solution .

We note that the convergence tested here is limited to the bubble filling factor, based on the radius of the outer shock and its fading time. We do not attempt to study convergence in other dynamical quantities such as the deposited momentum or the energy budget, as studied by Gentry et al. (2017).

7 Photo-ionized Bubbles

7.1 Introduction

Massive O/B stars, which precede the SNe, emit UV radiation that ionizes a Strömgren sphere around them, creating a bubble in which the SFR is suppressed, as in the SN bubble, but without totally evacuating the bubble interior. On galactic scales, supernovae clearly dominate the energy and momentum budget. Simulations show that photoionization alone cannot prevent the ionized ISM gas from catastrophically collapsing into molecular clouds and forming stars too rapidly, while the inclusion of SN feedback generates global turbulence in the ISM that properly suppresses the collapse to molecular clouds and maintains a diffuse atomic phase. The stellar-driven bubbles may dominate the destruction of molecular gas prior to the ignition of SNe in small star-forming clusters (Matzner, 2002; Krumholz & Matzner, 2009; Fall, Krumholz & Matzner, 2010), but in massive molecular clouds photo-ionization is expected to become ineffective once the escape velocity exceeds $5\,{\rm km}\,{\rm s}^{-1}$ (simulations by Dale, 2017). Under any conditions where photoionization dominates, one may apply in the photo-bubble case considerations involving the bubble filling factor analogous to the derivation of the KS relation by SN bubbles.

Consider an emission rate of ionizing photons ($h\nu>13.6$ eV) $Q=10^{49}\,Q_{49}$, and $T=10^{4}{\rm K}\,T_{4}$ for the equilibrium temperature within the ionized sphere777This temperature varies from $6\times 10^{3}$K to $1.5\times 10^{4}$K depending on the source and the medium. For a typical O star, with color temperature of $4\times 10^{4}$K, in a medium of solar metallicity and density of $100\,{\rm cm}^{-3}$, the equilibrium temperature is $\sim 7\times 10^{3}$K. In a very crude approximation, assuming a uniform density in a static bubble (e.g. Draine, 2011, Chapter 15), the Strömgren sphere radius $R_{\rm st}$ is determined by equilibrium between ionization and recombination rates, of the form $Q\propto R_{\rm st}^{3}n^{2}$, yielding $R_{\rm st}\propto n_{0}^{-2/3}$. If the relevant “fading” time is the fixed lifetime of the O star, the volume filling factor of the static Strömgren bubbles becomes

[TABLE]

where $S_{\rm st}=10^{-4}\,{\rm kpc}^{-3}\,{\rm yr}^{-1}{S_{\rm st}}_{-4}$ is the star-formation rate density of O stars. This would have implied a slope of $s=2$ in the KS relation.

However, the above static estimate does not account for the excess pressure created inside the Strömgren sphere, which makes the bubble expand beyond the static estimate. The pressure excess could range from a factor of two due to the doubling of the number of particles by the ionization (if the temperatures inside and outside the bubble are similar) to a factor of a few hundreds (if the outside medium is significantly cooler than the inside, as in a dense star-forming molecular cloud). Moreover, a spatially dependent radiation pressure from the star is also expected to affect the dynamics of the bubble. Therefore, one should appeal to a dynamical evaluation of the pressure inside and outside the bubble in order to properly estimate the final radius of the bubble, when the shell velocity has slowed down to the background speed of sound or turbulence velocity dispersion.

Following Krumholz & Matzner (2009), there is a characteristic radius (and time) below which the expansion is driven by radiation pressure and above which by gas pressure, approximated by (Krumholz & Matzner, 2009; Fall, Krumholz & Matzner, 2010)

[TABLE]

where the standard values for the other parameters have been assumed. This implies that the radiation pressure is effective only early in the central region while most of the evolution is dominated by the gas pressure, so we consider only the gas pressure below. Under the assumption that the shell expands slowly enough that the ionized gas has time to become uniform in density, which requires that sound waves be able to cross the ionized bubble, the self-similar expansion rate is

[TABLE]

[TABLE]

The radius and velocity of the shell at time $t=t_{6}\,{\rm Myr}$ are thus

[TABLE]

[TABLE]

This similarity solution is only valid for expansion speeds below the sound speed in the ionized gas $\sim 10\,{\rm km}\,{\rm s}^{-1}$). This assumption is satisfied or close to it for all real HII regions (Draine, 2011).

7.2 Photo-bubbles in two regimes

We consider a single star, or a short-duration cluster where $Q_{49}\sim\nu$, and where the lifetime of the ionizing source is $t_{*}=5\,{\rm Myr}\,\,t_{*,{\rm 5}}$. The fiducial value of $\sim 5\,{\rm Myr}$ is expected for massive stars in a cluster of $10^{6}M_{\odot}$ (Leitherer et al., 1999).

The stalling of the photo-bubbles can be evaluated in two different regimes of parameter space. In case A, the shell velocity slows down to the background sound speed, $v=c$, while the continuous ionization source is still alive, $t_{\rm st}<t_{*}$. In case B, $c<v(t_{*})<2c$, making the shell stall at $t_{*}$ (see below). In hypothetical case C, $v(t_{*})>2c$, the shell could have continued to expand in a snow-plow phase after $t_{*}$ until stalling with $v=c$ at $t_{\rm st}>t_{*}$. However, this is not likely to happen in HII regions, where the expansion speed is limited to $\leq 10\,{\rm km}\,{\rm s}^{-1}$ while the background sound speed (or turbulence velocity) is typically \lower 2.58333pt\hbox{;\buildrel>\over{\sim};}5\,{\rm km}\,{\rm s}^{-1}.

7.2.1 Case A: stall while the ionization source is alive

The shell expansion halts once the speed of the ionization front becomes equal to the external sound speed or turbulence velocity dispersion, while the ionization source is still alive. Based on eq. (77) and eq. (76), the stalling time and radius are

[TABLE]

[TABLE]

The volume filling factor of bubbles becomes888A full numerical solution including the radiation pressure term in eq. (74) indicates that this analytic estimate underestimates the filling factor by $\sim 30\%$ but the dependencies on the parameters are correct.

[TABLE]

For a constant $f_{\rm st}$ this implies a KS slope $s\simeq 2.6$ in case A.

The range of validity of case A, defined by $t_{\rm st}\leq t_{*}$, translates via eq. (78) to

[TABLE]

For a fiducial star cluster with a lifetime $t_{*,{\rm 5}}=1$, embedded in a molecular cloud of $n_{0}=100$ and $c_{1}=0.5$, this reads $Q_{49}\leq 300$, which refers to a low-mass star cluster of $M_{\rm s}\simeq 100\,Q_{49}\leq 3\times 10^{4}M_{\odot}$. Note that this upper-limit mass is very sensitive to the sound speed (or velocity dispersion).

7.2.2 Case B: stall when the ionization source dies

If the shell is still expanding with a velocity well above the background speed of sound once the ionization source shuts off, it will tend to enter a snow-plow expansion phase of $r\propto t^{2/7}$. In the transition from the expansion as $r\propto t^{4/7}$ to $r\propto t^{2/7}$, the velocity would drop by a factor of two. This implies that if $c<v(t=t_{*})<2c$, the expansion will stall right at $t_{*}$ with no further expansion.

Substituting the condition $v(t_{*})\geq 2c$ in eq. (77) gives

[TABLE]

For the fiducial star cluster above, this reads $Q_{49}\leq 3.2\times 10^{4}$, namely $M_{\rm s}\leq 3.2\times 10^{6}M_{\odot}$. Again, this upper-limit mass is very sensitive to the value of $c$.

Substituting for the stalling time $t_{\rm st}=t_{*}$ in eq. (76) gives

[TABLE]

[TABLE]

This implies a KS slope $s=0.86$ in case B.

7.2.3 Hypothetical Case C: stall after the ionization source dies

In the hypothetical case where $v(t_{*})>2c$ (possible only if $c<5\,{\rm km}\,{\rm s}^{-1}$), for a massive cluster where eq. (82) is invalid, when the ionizing photons shut off the shell could in principle enter a snow-plow phase, $r\propto t^{2/7}$. The stalling time and radius and the bubble filling would then be

[TABLE]

[TABLE]

[TABLE]

This implies a KS slope $s=1.6$ in this unlikely case C.

In summary, we expect for photo-bubbles alone (ignoring the following SNe) three regimes, depending on $Q_{49}$, namely on the star-cluster mass. For the fiducial values of the parameters ($t_{*,{\rm 5}}=1$, $n_{0}=100$ and $c_{1}=0.5$), cases A, B and C are expected to be valid for cluster masses $M_{\rm s}<3\times 10^{4}M_{\odot}$, $3\times 10^{4}<M_{\rm s}<3.2\times 10^{6}M_{\odot}$ and $M_{\rm s}>3.2\times 10^{6}M_{\odot}$, respectively, with case C being unrealistic (possibly valid only for ISM sound speeds lower than $5\,{\rm km}\,{\rm s}^{-1}$). The KS slopes are expected to be $s=2.6,0.86,1.6$ respectively. However, as we will see next (and as we saw in the simulations, §4), the photo bubbles are expected to be overwhelmed by the SN bubbles.

7.3 SNe in a photo-bubble

For a single O star followed by a single SN (or a short burst of stars), one can first envision a Strömgren bubble growing about the O star. During this period the SFR is suppressed within the bubble and the interior gas density is reduced due to the expansion of the shell. Then the SN turns on, generating a bubble that grows more-or-less following the standard evolution of a SN bubble but in a lower medium density. The reduced number density inside the photo-bubble, according to eq. 3 of Krumholz & Matzner (2009) and eq. (74), is

[TABLE]

assuming that the SN goes off at time $t_{*}$ after the birth of the O/B star. Substituting the density from eq. (88) in eq. (21) yields for the SN bubble filling factor

[TABLE]

To evaluate whether the SN bubble will grow bigger than the photo bubble, we compare the SN filling factor from eq. (89) to the photo-bubble filling factor in each of the three cases. Assuming $e_{51}=Q_{49}=\nu$, we obtain

[TABLE]

For the fiducial choice ($t_{*,{\rm 5}}=1$, $n_{0}=100$, $c_{1}=0.5$), we get $\nu\simeq 300$ and $3\times 10^{4}$ for the threshold values of $Q$ separating photo-bubble cases A from B and B from C respectively. With $\nu<300$ in case A we get $f_{0}/f_{\rm st}>104$, with $\nu>300$ in case B we get $f_{0}/f_{\rm st}>117$, and with $\nu>3\times 10^{4}$ in case C we get $f_{0}/f_{\rm st}>724$. This implies $f_{0}/f_{\rm st}\gg 1$ in all three cases, namely the SN bubbles are expected to be larger than the photo-bubbles and dominate the filling factor and therefore the KS relation.

Another factor that might have reduced the importance of the photo bubble is the accompanying stellar wind, which may take over the dynamics before the SN bubble takes over, as indicated in a 1D simulation of a $10^{6}M_{\odot}$ cluster by Gupta et al. (2016). They assumed an instantaneous burst of star formation and the corresponding sequence of SNe, where the radiation and mechanical luminosities of winds and SNe were computed using STARBURST99 (Leitherer et al., 1999). This cluster is in what we term zone B2.999With a SN-burst duration of $\sim 40\,{\rm Myr}$ and inter-SN interval $t_{\rm s}\sim 0.004\,{\rm Myr}$, and with a cooling time for an individual SN at $n_{0}\sim 100$ of $t_{\rm c}\sim 0.005\,{\rm Myr}$ and fading time $t_{\rm f}\sim 25\,{\rm Myr}$ for the super-bubble, we have $t_{\rm s}<t_{\rm c}$ and $t_{\rm f}<t_{\rm b}$. Almost independent of background density, the radiation force turned out to be important only during the first Myr, after which the stellar wind takes over in this simulation, and the photo-ionized bubble collapses to the wind-driven thin shell. Eventually, after $\sim 3\,{\rm Myr}$, the dynamics of the super-bubble is driven by the SNe. One should note, however, that the importance of the wind in the intermediate stage is controversial, as it depends criticality on how efficiently it is trapped in the HII region, versus how much of it leaks out (Krumholz, McKee & Bland-Hawthorn, 2018). In simulations, this depends sensitively on the initial conditions, and in 1D simulations it depends on the subgrid model used to represent wind leakage. Observations suggest efficient leakage, which makes winds sub-dominant compared to ionized gas and radiation (Lopez et al., 2011).

We conclude that for a short burst of clustered stars and SNe the SN filling factor is expected to dominate the bubble filling factor and thus determine the KS relation. This is confirmed in the isolated-galaxy simulations presented in §4, where Fig. 8, Fig. 9 and Fig. 10 demonstrate that the SN-feedback overwhelms the HII feedback as the mechanism that produces the constant hot filling factor and provides a tight global KS relation with a global slope $s\simeq 1.5$. It should be borne in mind, however, that this encouraging confirmation of our theoretical considerations is naturally limited to the local SFR recipes and feedback processes as implemented in these specific simulations.

8 Discussion

Here we address some of our key assumptions, discuss potential caveats, and comment on the observed KS relation.

8.1 Negative or positive feedback?

Among our simplifying assumptions, we adopted here the notion that SN feedback (and any other stellar feedback) is negative, suppressing the SFR in the SN bubbles. This is motivated by the robust fact that when SN feedback is incorporated in simulations, the timescale for star formation typically slows down from the free-fall time ($\epsilon_{\rm ff}\sim 1$) to a timescale larger by two orders of magnitude ($\epsilon_{\rm ff}\sim 0.01$). We thus ignore here the counter possibility that stellar or SN feedbacks may actually be positive, triggering star formation in the bubble shells. This has been addressed in a limited way by theory (Elmegreen & Lada, 1977) and observations (Samal et al., 2014; Egorov et al., 2017), the latter showing indications for high SFR near feedback-driven super-bubble structures from O-type stars. However, based on simulations, caution has been called for when interpreting these observations in terms of “triggering” (Dale, Haworth & Bressert, 2015). It is obvious that SFR is suppressed in the interior of a SN bubble, which is hot and dilute. However, after a short while, most of the bubble gas has been swept into a dense and relatively cold shell, raising the question of why the SFR is suppressed there. The common wisdom is that the suppression is largely by mechanical effects, such as ram pressure associated with the high velocity of the shell with respect to the ISM, or the high turbulence within the very perturbed and non-uniform shell and the associated strong shear. This was partly demonstrated, e.g., in high-resolution simulations of SNe in a single molecular cloud (Rogers & Pittard, 2013). While the issue of negative versus positive feedback is still an open issue which deserves further study, as it is for AGN feedback (Silk, 2013; Bieri et al., 2015), our analysis was based on the assumption that the SFR is suppressed in all the ISM gas that has been swept by the SN bubble. The swept gas mass filling factor can thus be approximated by the volume filling factor of the hot bubbles.

8.2 Molecular gas

In our model we have assumed that the portion of the ISM that is not filled by SN bubbles is mostly star-forming molecular hydrogen. We have not included the possibility that a significant portion of the ISM might consist of non-star-forming, warm HI. As applied to weakly star-forming regions such as dwarf galaxies and the outer disks of modern spiral galaxies, this is in fact not a good assumption. In these regions warm HI dominates the neutral portion of the ISM by both mass and volume, and star formation is likely regulated at least in part by the thermal and chemical transitions between the warm, non-self-gravitating, non-star-forming and cold, bound, star-forming phases (Krumholz, McKee & Tumlinson, 2009b; Ostriker, McKee & Leroy, 2010; Krumholz, 2013; Kim, Ostriker & Kim, 2013; Forbes et al., 2016) For such galaxies the observed KS relation is significantly more complex than the simple power law relation and constant $\epsilon_{\rm ff}$ that describes the molecular phase; the star formation rate in atomic gas shows significant dependence on third quantities such as metallicity and stellar density (Bolatto et al., 2011; Shi et al., 2014; Jameson et al., 2016). The simple model we present here is not intended to apply to galaxies with a significant warm atomic component.

Fortunately, this is a modest limitation when it comes to high-redshift galaxies. Both observations (Bigiel et al., 2008; Lee et al., 2012; Wong et al., 2013) and theory (Krumholz, McKee & Tumlinson, 2009a; Sternberg et al., 2014) suggest that galaxies transition from atomic-dominated to molecule-dominated at a surface density $\approx 10(Z/Z_{\odot})^{-1}\,M_{\odot}\,{\rm pc}^{-2}$, where $Z$ is the metallicity. The VELA simulations that we use to test our model are all well above this threshold (c.f. Fig. 2), as are essentially all observed high-redshift galaxies (Genzel et al., 2010; Daddi et al., 2010a; Tacconi et al., 2013). Thus our model applies to the majority of observable star-forming systems beyond the local Universe.

For low-redshift galaxies the situation is somewhat more complex. Local dwarf galaxies have surface densities below the critical $\approx\!10(Z/Z_{\odot})^{-1}M_{\odot}\,{\rm pc}^{-2}$ that marks the HI to H2 transition throughout their area, and thus are almost entirely HI-dominated. Local spirals, on the other hand, tend to cross the threshold from mostly HI to mostly H2 at radii $\approx\!0.5\,r_{25}$ (Leroy et al., 2008; Schruba et al., 2011), a radius that is relatively close to the stellar scale length. Star formation follows H2, so this is also the radius within which most star formation is concentrated. Our model is therefore reasonably applicable to the central parts of modern spirals, where a substantial fraction of their stars form.

8.3 Potential effect of gravity

Our analytic modeling ignores the gravitational force acting on the expanding bubbles, which may, in principle, affect the expansion and even cause a re-collapse. This is discussed using sophisticated 1D models by Rahner et al. (2018, 2019). However, their analysis of the re-collapse effect is mostly concerned with non-SN feedback. For SNe, this is is unlikely to be important because of the following argument. A single SN, even if any enhancement from clustering is ignored, reaches a terminal momentum of $\sim 3\times 10^{5}M_{\odot}\,{\rm km}\,{\rm s}^{-1}$. Averaging over an IMF that makes 1 SN per $100M_{\odot}$, this gives $3\times 10^{3}M_{\odot}\,{\rm km}\,{\rm s}^{-1}$ per $1M_{\odot}$ of stars formed. If $1\%$ of the cloud mass is turned into stars, the momentum budget per mass of gas in the GMC is reduced to only $30\,{\rm km}\,{\rm s}^{-1}$ This is to be compared to the typical escape speeds of GMCs, which are closer to $5\,{\rm km}\,{\rm s}^{-1}$. This indicates that self-gravity is not a relevant consideration for SN remnant evolution. It might be relevant for other feedback mechanisms, which are working with much lower momentum budgets.

8.4 Super-bubble blow-out

The analysis so far assumed that the bubbles are confined to the galactic disc. This is probably a valid assumption for unclustered SNe, where the fading radius is on the order of a few tens of parsecs (eq. 19), significantly smaller than the typical disc height. However, super-bubbles of multiple SNe in massive clusters (high $\nu$) could have a fading radius that exceeds the disc height (eqs. 45, 60, 55). In this case, the bubbles will blow out from the disc before completing their expansion. Cold gas from the shells is expected to be ejected with velocities comparable to the escape velocity from the disc, while the over-pressured hot gas from the bubble interior is expected to escape at higher velocities. This may change the considerations concerning the self-regulation to a constant bubble filling factor, which should now consider the balance between outflows and inflows, including recycling. Self-regulation could in principle be achieved in this case as well, as a high SFR would be associated with strong outflows, which remove (mostly cold) gas and thus decrease the SFR. To maintain the SFR, the disc has to be fueled with cold gas, by cosmological inflow or by return of the outflowing cold gas. Alternatively, it is possible that the bubble filling factor is self-regulated by bubble blow-out without involving the SFR and SN feedback at all, namely once the filling factor is larger than $\sim 0.5$, it generates outflows from the disc that keeps $f$ near $0.5$. We defer the analysis of these cases of major outflows and inflows to future work.

If the outflows are non-negligible but they do not change the overall picture where the bubbles are largely confined to the disc, one may repeat the analysis with the fading radius replaced by the half-height of the disc, $h$, assuming that the bubble stops growing within the disc plane once it blows out of the disc in the vertical direction. In this case the relevant bubble volume is $V_{\rm sn}\sim h^{3}$. The disc height is determined by the vertical balance between gravity and pressure gradient, $dP/dz\propto G\rho\Sigma$, where $\Sigma$ is the surface density in the disc. With $dP/dz\propto\rho c^{2}/h$, this gives $h\propto cn^{-1/2}$, or

[TABLE]

If we crudely assume that the time for the bubble to reach the disc edge scales like the single-SN fading time, $t(h)\propto n^{-0.37}$, the bubble filling factor would be

[TABLE]

With a constant $f$ we therefore expect a KS relation $\rho_{\rm sfr}\propto n^{s}$ with $s$ in the range $1.5-2.0$. Alternatively, for clustered SNe, one can use the shell expansion rate in the continuous source phase, $R_{\rm s}\propto t^{3/5}\,n^{-1/5}$ (eq. 30 or eq. 52). Assuming as above $R_{\rm s}=h$ and $h\propto cn^{-1/2}$, one obtains $t(h)\propto c^{1.67}\,n^{-0.5}$, and the bubble filling factor becomes

[TABLE]

namely $s=2$ in the KS relation.

The importance of outflows can be evaluated by comparing the timescale for outflows with the timescale for self-regulation of the bubble filling factor. For the latter we take as a reference the bubble fading time, which is $t_{\rm f}\sim 1\,{\rm Myr}$ for single SNe, $t_{\rm f}\sim 10\,{\rm Myr}$ for SNe in low-mass clusters, and $t_{\rm f}\sim 100\,{\rm Myr}$ for massive clusters. Consider Milky-Way-like galaxies as an example, where the outflow rate, if comparable to the SFR, is \dot{M}\lower 2.58333pt\hbox{;\buildrel<\over{\sim};}10M_{\odot}\,{\rm yr}^{-1}, and the gas mass is $M\sim 10^{10}M_{\odot}$. The timescale for gas depletion by outflow is thus t_{\rm out}=M/\dot{M}\lower 2.58333pt\hbox{;\buildrel>\over{\sim};}1\,{\rm Gyr} (or even longer for the large-scale outflow, given that it contributes only a fraction of the total outflow). This is much longer than the SN fading time, indicating that outflows from Milky-Way-type discs are not expected to have a strong effect on the filling-factor considerations. Our isolated-galaxy simulations indeed show similar outflow timescales. For example, the galaxy with 15% gas fraction (mimicking $z=0$) has $t_{\rm out}\simeq 1.6\,{\rm Gyr}$ (and $\dot{M}_{\rm out}\simeq 0.7\,{\rm SFR}$), while the galaxy with 50% gas ($z\sim 2$) has $t_{\rm out}\simeq 300\,{\rm Myr}$.

In the VELA cosmological simulations at $z\sim 1$, we measure the typical outflow timescales to be $t_{\rm out}\sim 100\,{\rm Myr}$. This is consistent with the typical SFR and cosmological accretion rate being higher by an order of magnitude at $z\sim 1$ compared to $z=0$. This is still longer than the fading time for single SNe and for SNe in low-mass clusters, but it is comparable to $t_{\rm f}$ for massive star clusters, indicating that SN blow-outs should be considered in the case of massive star-forming clusters.

One may ask whether accretion onto the disc may affect the considerations. The typical cosmological accretion timescale (in the EdS regime, $z>1$ say) is (Dekel et al., 2013)

[TABLE]

This is larger than $1\,{\rm Gyr}$ for $z<3$, namely it is much longer than the SN fading time, indicating that the effect of cosmological accretion is not directly relevant for our purposes. On the other hand, the effect of recycled gas, returning after it outflew from the disc, may be more important, as indicated by comparing theoretical and observational specific SFR at $z=1-4$ (e.g. Dekel & Mandelker, 2014). Assuming an initial vertical outflow velocity on the order of the circular velocity, the return time would be comparable to the disc dynamical time, typically on the order of $100\,{\rm Myr}$. This indicates that the recycling should affect the disc gas budget in the case of SNe in massive star clusters.

8.5 On the observed KS relation

While our target seemed to be the slope $s\simeq 1.5$ of the macroscopic relation, we note that the observed value of the slope may deviate from $1.5$, in the range $s=1-2$. For example, in patches of $\sim 1\,{\rm kpc}^{2}$ in nearby galaxies, the typical obtained slopes are $s\simeq 1$ (Leroy et al., 2013, 2017; Utomo et al., 2018). On the other hand, referring to whole-galaxy averages, and including high-redshift and starburst galaxies, the estimates range from $s\sim 1$ (Genzel et al., 2011) to $s\sim 2$ (Faucher-Giguère, Quataert & Hopkins, 2013), partly depending on the assumed $\alpha_{\rm CO}$, the CO-to-H2 conversion factor. A more complete data compilation consistently yields slopes in the range $s=1-2$ (Krumholz et al., 2017).

As another example, Tacconi et al. (2018), who evaluated gas masses $M_{\rm g}$ for hundreds of galaxies in the redshift range $z=0-4$, obtained a gradient of the depletion time across the Main Sequence of star-forming galaxies, at given stellar mass and redshift,

[TABLE]

where $\dot{M}_{\rm sf}$ is the star-formation rate. The best-fit slope from the data is $\tau=0.44$ but the minimum uncertainty of $\pm 0.16$ dex, largely reflecting systematics in the different analyses, allows values of $\tau$ from below $0.3$ to above $0.6$. The deduced KS relation is $\dot{M}_{\rm sf}\propto M_{\rm g}^{s}$ with $s=1/(1-\tau)$, namely, the best-fit KS slope is $s=1.79$, but values from below $1.4$ to above $2.0$ can be accommodated. We note that this data may mix relatively relaxed discs and (merger-induced?) starbursts.

9 Conclusion

We suggest that the relation between the macroscopic densities of SFR and gas mass, the KS relation, may naturally arise from the way supernova bubbles evolve, independent of the details of the microscopic star-formation law, and with no explicit dependence on gravity.

The key idea is that the filling factor $f$ of the SN bubbles in which SFR is suppressed is self-regulated into a constant value of order one half; a larger (smaller) filling factor causes a slowdown (speedup) of the SFR and hence the SN rate, which reduces (increases) the filling factor. This has been demonstrated using both zoom-in cosmological simulations and isolated-galaxy simulations, using different codes, subgrid recipes and initial/boundary conditions.

Given the bubble fading radius and time, $R_{\rm f}(n)$ and $t_{\rm f}(n)$, as a function of hydrogen number density in the ISM, $n$, the filling factor of gas where SFR is suppressed can be expressed in terms of the SN rate density, $S$, as $f\propto S\,n^{-s}$. With $f$ fixed at a constant value, this implies for the macroscopic SFR density $\rho_{\rm sfr}\propto S\propto n^{s}$.

An analytic toy model based on the standard evolution of SN bubbles, assumed to be spread at random in space, predicts a slope of $s=1.48$, suspiciously close to the observed value. We generalized the toy-model analysis to a sequence of SNe coincidental in space and of different patterns in time, exploring different regimes in parameter space, using analytic modeling and spherical simulations. While the analytic predictions for the slope range from $s=1$ to $s=2$ in different cases, the slopes in the simulations tend to be closer to $s\simeq 1.5$ in most cases.

When expressed in terms of the free-fall time in macroscopic volumes, $\rho_{\rm sfr}=\epsilon_{\rm ff}\rho_{\rm g}/t_{\rm ff}$, the model predicts that the efficiency factor $\epsilon_{\rm ff}$ is independent of $n$, and is of order $0.01-0.02$, in the ball park of the observed values, both on microscopic and macroscopic scales.

An analogous analytic toy model for suppressed SFR in photo-ionized bubbles around the massive stars that precede the SNe indicates slopes in a broader range depending on the mass of the star-forming cluster. However, the ionized spheres collapse to the wind-driven super-bubbles after the first Gigayear, and the super-baubles are eventually dominated by the SN energy for most of their lifetimes, indicating that the SN super-bubbles dominate the filling factor and determine the KS relation.

Our zoom-in cosmological simulations as well as our isolated-disc simulations convincingly demonstrate the self-regulation to a constant hot filling factor in realistic galactic discs, and the associated generation of a KS relation. The isolated-disc simulations, run with different local SFR recipes, demonstrate that the self-regulation to a constant hot filling factor is indeed insensitive to the local SFR law, and so is the resultant KS relation. The latter is consistent with the simulated results of Hopkins, Quataert & Murray (2011) and Hopkins, Narayanan & Murray (2013) that the global SFR is robust to variations in the local SFR recipes over a broad range of local recipes. Our simulations also indicate that the self-regulation is due to feedback, and that SN feedback dominates over photo-ionization feedback.

Our analytic modeling of the bubbles and their filling factor is clearly a crude over-simplification in many ways. For example, we ignore the complex spatial structure of star-forming clouds, the inhomogeneous nature of the ISM and the effects of neighboring bubbles running into each other. We address only in a crude way the interplay between the effects of pre-SN massive stars, through winds, photo-ionization and radiation pressure, and the evolution of SN bubbles. Among other shortcomings, we do not deal with the formation of molecular gas, ignore the ejection of gas from the galactic disc, and neglect the possible positive aspects of the feedback. For the global evolution of the discs, we do not consider the accretion onto the disc, the self-regulation of disc instability, the associated instability-driven inflow within the disc, the generation of turbulence by the above and by feedback, and the vertical balance between turbulence and gravity. Despite all these shortcomings, the simple concept of self-regulation to a constant filling factor for the gas where star formation is suppressed, combined with the evolution of SN bubbles based on common physics, gives rise in a trivial way to a universal relation between the macroscopic SFR and gas content in galaxies or in patches within galaxies, similar to the observed KS relation. Given the severe shortcomings of our model one may argue that it is too simple to be true. On the other hand, its simplicity and robustness may indicate that it captures the essence of the origin of the KS relation.

The confirmed validity of the self-regulation into a constant filling factor in the simulations, cosmological and isolated, and the demonstrated robustness of the macroscopic relation to the local star-formation recipe, provide encouraging evidence for the validity of this basic scenario. It argues that supernova feedback is the key factor in the origin of the KS relation.

Future work may proceed in several routes. First, one should improve the idealized analytic modeling for better agreement between the analytic predictions and the spherical simulations. Then, one may attempt to generalize the modeling to address some of the effects that were ignored in the current simplified analysis. In particular, one should consider the alternative possibility of self-regulation by outflows and inflows. Finally, one should analyze simulated galaxies in greater detail, e.g., to follow the evolution of super-bubbles and the way the self-regulation to a constant filling factor is achieved.

Acknowledgments

This work was inspired by an interaction with Jerry Ostriker. We are grateful for stimulating discussions with Siddhartha Gupta, Miao Li, Chris McKee, Brant Robertson and Amiel Sternberg. This work was partly supported by the grants France-Israel PICS, Germany-Israel GIF I-1341-303.7/2016, Germany-Israel DIP STE1869/2-1 GE625/17-1, I-CORE Program of the PBC/ISF 1829/12, ISF 857/14, US-Israel BSF 2014-273, and NSF AST-1405962. The cosmological VELA simulations were performed at the National Energy Research Scientific Computing Center (NERSC) at Lawrence Berkeley National Laboratory, and at NASA Advanced Supercomputing (NAS) at NASA Ames Research Center. Development and analysis have been performed in the astro cluster at HU. The isolated-galaxy simulations were performed using the HPC resources of CINES and TGCC under the allocations A0030402192 and A0050402192 made by GENCI.

Appendix A The VELA Cosmological Simulations

The VELA suite consists of hydro-cosmological simulations of 35 moderately massive galaxies. Full details are presented in Ceverino et al. (2014); Zolotov et al. (2015). This suite has been used to study central issues in the evolution of galaxies at high redshifts, including compaction to blue nuggets and the trigger of quenching (Zolotov et al., 2015; Tacchella et al., 2016b, a), evolution of global shape (Ceverino et al., 2015; Tomassetti et al., 2016), violent disc instability (Mandelker et al., 2014, 2017), OVI in the CGM (Roca-Fàbrega et al., 2018), and galaxy size and angular momentum (Jiang et al., 2018). Additional analysis of the same suite of simulations are discussed in Moody et al. (2014); Snyder et al. (2015). In this appendix we give an overview of the key aspects of the simulations and their limitations.

A.1 The Cosmological Simulations

The VELA simulations make use of the Adaptive Refinement Tree (ART) code (Kravtsov, Klypin & Khokhlov, 1997; Kravtsov, 2003; Ceverino & Klypin, 2009), which accurately follows the evolution of a gravitating N-body system and the Eulerian gas dynamics using an adaptive mesh refinement approach. The adaptive mesh refinement maximum resolution is $17-35\,{\rm pc}$ at all times, which is achieved at densities of $\sim 10^{-4}-10^{3}\,{\rm cm}^{-3}$. Beside gravity and hydrodynamics, the code incorporates physical process relevant for galaxy formation such as gas cooling by atomic hydrogen and helium, metal and molecular hydrogen cooling, photoionization heating by the UV background with partial self-shielding, star formation, stellar mass loss, metal enrichment of the ISM and stellar feedback. Supernovae and stellar winds are implemented by local injection of thermal energy as described in Ceverino & Klypin (2009); Ceverino, Dekel & Bournaud (2010) and Ceverino et al. (2012). Radiation-pressure stellar feedback is implemented at a moderate level, following Dekel et al. (2013), as described in Ceverino et al. (2014).

Cooling and heating rates are tabulated for a given gas density, temperature, metallicity and UV background based on the CLOUDY code (Ferland et al., 1998), assuming a slab of thickness $1\,{\rm kpc}$. A uniform UV background based on the redshift-dependent Haardt & Madau (1996) model is assumed, except at gas densities higher than $0.1\,{\rm cm}^{-3}$, where a substantially suppressed UV background is used ($5.9\times 10^{6}\,{\rm erg}\,{\rm s}^{-1}\,{\rm cm}^{-2}\,{\rm Hz}^{-1}$) in order to mimic the partial self-shielding of dense gas, allowing dense gas to cool down to temperatures of $\sim 300$K. The assumed equation of state is that of an ideal mono-atomic gas. Artificial fragmentation on the cell size is prevented by introducing a pressure floor, which ensures that the Jeans scale is resolved by at least 7 cells (see Ceverino, Dekel & Bournaud, 2010).

Star particles form in timesteps of $5\,{\rm Myr}$ in cells where the gas density exceeds a threshold of $1~{}\,{\rm cm}^{-3}$ and the temperatures is below $10^{4}$K. Most stars ($>90\%$) end up forming at temperatures well below $10^{3}$K, and more than half of the stars form near $300$K in cells where the gas density is higher than $10~{}\,{\rm cm}^{-3}$. The code implements a stochastic star-formation where a star particle with a mass of $42\%$ of the gas mass forms with a probability $P=(\rho_{\rm g}/10^{3}\,{\rm cm}^{-3})^{1/2}$ but not higher than $0.2$. This corresponds to a local SFR that crudely mimics $\rho_{\rm sfr}\epsilon_{\rm ff}\rho_{\rm g}/t_{\rm ff}$ with $\epsilon_{\rm ff}\sim 0.02$. A stellar initial mass function of Chabrier (2003) is assumed.

Thermal feedback that mimics the energy release from stellar winds and supernova explosions s incorporated as a constant heating rate over the $40~{}\,{\rm Myr}$ following star formation. A velocity kick of $\sim 10\,{\rm km}\,{\rm s}^{-1}$ is applied to $30~{}\%$ of the newly formed stellar particles – this enables SN explosions in lower density regions where the cooling may not overcome the heating without implementing an artificial shutdown of cooling (Ceverino & Klypin, 2009). The code also incorporates the later effects of Type Ia supernova and stellar mass loss, and it follows the metal enrichment of the ISM.

Radiation pressure is incorporated through the addition of a non-thermal pressure term to the total gas pressure in regions where ionizing photons from massive stars are produced and may be trapped. This ionizing radiation injects momentum in the cells neighbouring massive star particles younger than $5\,{\rm Myr}$, and whose column density exceeds $10^{21}\,{\rm cm}^{-2}$, isotropically pressurizing the star-forming regions (see more details in Agertz et al., 2013; Ceverino et al., 2014).

The initial conditions for the simulations are based on dark-matter haloes that were drawn from dissipationless N-body simulations at lower resolution in cosmological boxes of $15-60\,{\rm Mpc}$. The $\Lambda$CDM cosmological model was assumed with the WMAP5 values of the cosmological parameters, $\Omega_{\rm m}=0.27$, $\Omega_{\Lambda}=0.73$, $\Omega_{\rm b}=0.045$, $h=0.7$ and $\sigma_{8}=0.82$ (Komatsu et al., 2009). Each halo was selected to have a given virial mass at $z=1$ and no ongoing major merger at $z=1$. This latter criterion eliminated less than $10~{}\%$ of the haloes, those that tend to be in a dense, proto-cluster environment at $z\sim 1$. The virial masses at $z=1$ were chosen to be in the range $M_{\rm v}=2\times 10^{11}-2\times 10^{12}~{}M_{\odot}$, about a median of $4.6\times 10^{11}~{}M_{\odot}$. If left in isolation, the median mass at $z=0$ was intended to be $\sim 10^{12}~{}M_{\odot}$.

The VELA cosmological simulations are state-of-the-art in terms of high-resolution adaptive mesh refinement hydrodynamics and the treatment of key physical processes at the subgrid level. In particular, they trace the cosmological streams that feed galaxies at high redshift, including mergers and smooth flows, and they resolve the violent disc instability that governs high-$z$ disc evolution and bulge formation (Ceverino, Dekel & Bournaud, 2010; Ceverino et al., 2012, 2015; Mandelker et al., 2014). Like in other simulations, the treatments of star formation and feedback processes are rather simplified. The code may assume a realistic SFR efficiency per free fall time on the grid scale but it does not follow in detail the formation of molecules and the effect of metallicity on SFR. The feedback is treated in a crude way, where the resolution does not allow the capture of the Sedov-Taylor phase of supernova bubbles. The radiative stellar feedback assumed no infrared trapping, in the spirit of low trapping advocated by Dekel & Krumholz (2013) based on Krumholz & Thompson (2013), which makes the radiative feedback weaker than in other simulations that assume more significant trapping (Murray, Quataert & Thompson, 2010; Hopkins, Quataert & Murray, 2012). Finally, AGN feedback, and feedback associated with cosmic rays and magnetic fields, are not yet implemented. Nevertheless, as shown in Ceverino et al. (2014), the star formation rates, gas fractions, and stellar-to-halo mass ratio are all in the ballpark of the estimates deduced from observations.

A.2 The Galaxy Sample and Measurements

The virial and stellar properties of the galaxies are listed in Table 3. The virial mass $M_{\rm v}$ is the total mass within a sphere of radius $R_{\rm v}$ that encompasses an overdensity of $\Delta(z)=[18\pi^{2}-82\Omega_{\Lambda}(z)-39\Omega_{\Lambda}(z)^{2}]/\Omega_{\rm m}(z)$, where $\Omega_{\Lambda}(z)$ and $\Omega_{\rm m}(z)$ are the cosmological parameters at $z$ (Bryan & Norman, 1998; Dekel & Birnboim, 2006). The stellar mass $M_{\rm s}$ is measured within a radius of $0.2R_{\rm v}$.

We start the analysis at the cosmological time corresponding to expansion factor $a=0.125$ (redshift $z=7$). As can be seen in Table 3, most galaxies reach $a=0.50$ ($z=1$). Each galaxy is analyzed at output times separated by a constant interval in $a$, $\Delta a=0.01$, corresponding at $z=2$ to $\sim 100~{}\,{\rm Myr}$ (roughly half an orbital time at the disc edge). The sample consists of totally $\sim 1000$ snapshots in the redshift range $z=6-1$ from 35 galaxies that at $z=2$ span the stellar mass range $(0.2-6.4)\times 10^{11}M_{\odot}$. The half-mass sizes $R_{\rm e}$ are determined from the $M_{\rm s}$ that are measured within a radius of $0.2R_{\rm v}$ and they range $R_{\rm e}\simeq 0.4-3.2\,{\rm kpc}$ at $z=2$.

The SFR for a simulated galaxy is obtained by ${\rm SFR}=\langle M_{\star}(t_{\rm age}<t_{\rm max})/t_{\rm max}\rangle_{t_{\rm max}}$, where $M_{\rm s}(t_{\rm age}<t_{\rm max})$ is the mass at birth in stars younger than $t_{\rm max}$. The average $\langle\cdot\rangle_{t_{\rm max}}$ is obtained by averaging over all $t_{\rm max}$ in the interval $[40,80]\,{\rm Myr}$ in steps of $0.2\,{\rm Myr}$. The $t_{\rm max}$ in this range are long enough to ensure good statistics. The SFR ranges from $\sim 1$ to $33M_{\odot}\,{\rm yr}^{-1}$ at $z\sim 2$.

The instantaneous mass of each star particle is derived from its initial mass at birth and its age using a fitting formula for the mass loss from the stellar population represented by the star particle, according to which 10%, 20% and 30% of the mass is lost after 30 Myr, 260 Myr , and 2 Gyr from birth, respectively. We consistently use here the instantaneous stellar mass, $M_{\rm s}$, and define the specific SFR by ${\rm sSFR}={\rm SFR}/M_{\rm s}$.

The determination of the centre of the galaxy is outlined in detail in Appendix B of Mandelker et al. (2014). Briefly, starting form the most bound star, the centre is refined iteratively by calculating the centre of mass of stellar particles in spheres of decreasing radii, updating the centre and decreasing the radius at each iteration. We begin with an initial radius of 600 pc, and decrease the radius by a factor of 1.1 at each iteration. The iteration terminates when the radius reaches 130 pc or when the number of stellar particles in the sphere drops below 20.

The disc plane and dimensions are determined iteratively, as detailed in Mandelker et al. (2014). The disc axis is defined by the angular momentum of cold gas ($T<1.5\times 10^{4}$K), which on average accounts for $\sim 97\%$ of the total gas mass in the disc. The radius $R_{\rm d}$ is chosen to contain $85\%$ of the cold gas mass in the galactic mid-plane out to $0.15R_{\rm v}$, and the half-height $H_{\rm d}$ is defined to encompass $85\%$ of the cold gas mass in a thick cylinder where both the radius and half-height equal $R_{\rm d}$.

Relevant global properties of the VELA 3 galaxies at $z=2$ are listed in Table 3 and explained in the caption. It includes the global masses and sizes of the different components, and the quantities within the thin discs analyzed that are relevant for the KS relation and the hot filling factor.