Enhanced momentum feedback from clustered supernovae

Gentry, Eric S.; Krumholz, Mark R.; Dekel, Avishai; Madau, Piero

doi:10.1093/mnras/stw2746

Abstract

Young stars typically form in star clusters, so the supernovae (SNe) they produce are clustered in space and time. This clustering of SNe may alter the momentum per SN deposited in the interstellar medium (ISM) by affecting the local ISM density, which in turn affects the cooling rate. We study the effect of multiple SNe using idealized 1D hydrodynamic simulations which explore a large parameter space of the number of SNe, and the background gas density and metallicity. The results are provided as a table and an analytic fitting formula. We find that for clusters with up to ∼100 SNe, the asymptotic momentum scales superlinearly with the number of SNe, resulting in a momentum per SN which can be an order of magnitude larger than for a single SN, with a maximum efficiency for clusters with 10–100 SNe. We argue that additional physical processes not included in our simulations – self-gravity, breakout from a galactic disc, and galactic shear – can slightly reduce the momentum enhancement from clustering, but the average momentum per SN still remains a factor of 4 larger than the isolated SN value when averaged over a realistic cluster mass function for a star-forming galaxy. We conclude with a discussion of the possible role of mixing between hot and cold gas, induced by multidimensional instabilities or pre-existing density variations, as a limiting factor in the build-up of momentum by clustered SNe, and suggest future numerical experiments to explore these effects.

hydrodynamics, ISM: bubbles, ISM: supernova remnants

1 INTRODUCTION

Supernovae (SNe) play a key role in regulating star formation at galactic scales. SN energy, if retained, can disrupt molecular clouds and small galaxies (Dekel & Silk 1986). Even if significant energy is lost to radiative cooling, SNe inject momentum which cannot be radiated away, which drives turbulence, the dominant form of dynamical pressure support in galactic discs (Jenkins & Tripp 2011; Kim, Kim & Ostriker 2011). This turbulent support both prevents the collapse of star-forming regions (locally limiting star formation; Ostriker & Shetty 2011; Faucher-Giguère, Quataert & Hopkins 2013) and drives galactic winds (globally limiting star formation; Murray, Quataert & Thompson 2005; Hopkins, Quataert & Murray 2012; Creasey, Theuns & Bower 2013; Dekel & Krumholz 2013; Thompson & Krumholz 2016).

Unfortunately, the processes controlling supernova remnant (SNR) evolution operate at smaller scales than what can typically be resolved by galaxy or cosmological simulations. In particular, the dense shells of SNRs rapidly radiate away most of the SN energy, leaving a cold dense shell and a hot diffuse interior. If a simulation cannot resolve these two zones, then it cannot realistically evolve the SNR, resulting in problems such as overcooling (Katz 1992). To counteract this, some authors have prescribed turning off cooling for young, unresolved SNRs (Gerritsen 1997; Stinson et al. 2006), while others have proposed models which mimic the otherwise unresolved multiphase nature of the interstellar medium (ISM; Keller et al. 2014). These methods have their strengths, but the most direct way to incorporate the relevant physics is multiscale modelling: evolve a number of SNRs in a representative set of environments using simulations with high enough resolution to resolve the relevant processes, and use those results in large, low-resolution simulations.

Early attempts at using high-resolution simulations to create subgrid SN feedback models focused on producing energy-driven models (Thornton et al. 1998), but recently there has been increased interest in momentum-driven models, both by those using high-resolution simulations to study SNRs directly (Iffrig & Hennebelle 2015; Kim & Ostriker 2015; Martizzi, Faucher-Giguère & Quataert 2015; Walch & Naab 2015) and by those who use such models in lower resolution simulations (Hopkins, Quataert & Murray 2011; Shetty & Ostriker 2012; Hennebelle & Iffrig 2014; Goldbaum, Krumholz & Forbes 2016). This change in emphasis has been driven by the realization that, while the energy content of SNRs is important for producing hot galactic winds which are observable in X-rays, the momentum budget is more important when it comes to SNRs regulating star formation and possibly ejecting cool gas from galaxies (Dekel & Krumholz 2013).

At early times, before radiative losses are important, an SNR is in the Sedov stage, during which the energy is approximately conserved and the radial momentum is increasing. Once radiative losses become significant, it enters a pressure-driven snowplough phase, during which the energy is decreasing and the momentum is still increasing. As the bubble expands and cools (adiabatically and radiatively), its pressure will eventually decrease to the ISM pressure, at which point it becomes a momentum-driven snowplough. Asymptotically, in the idealized case of a spherical SNR expanding into a uniform, cold medium, this results in zero energy being added to the ISM, but a non-zero and finite amount of momentum being added. The goal of high-resolution simulations is to follow all of these phases, and identify the asymptotic momentum as a function of the properties of the driving stars and the large-scale environment, making this value available for use in larger scale models.

A number of authors have performed systematic parameter studies of the expansion of an SNR from a single SN in spherical symmetry (Chevalier 1974; Cioffi, McKee & Bertschinger 1988; Thornton et al. 1998). The most complete of these studies, that of Thornton et al. (1998), spanned metallicities from 10⁻³ to 10^0.5 times solar and ambient densities from 0.1 to 10³ H atoms cm⁻³. More recently, there have been a number of 3D simulations which allowed the study of SNR evolution within a more realistic, non-spherically symmetric background. Martizzi et al. (2015), Kim & Ostriker (2015), and Walch & Naab (2015) all found that inhomogeneities present prior to the first SN explosion – such as those expected due to a multiphase structure of the ISM or ionized bubbles created by pre-SN radiation – did not change the final momentum by more than 60 per cent. A more interesting effect was found by considering the inhomogeneities which result from bubbles of previous SNRs. Kim & Ostriker (2015) found that a series of clustered SNe can decrease the momentum per SN, in some cases by almost a factor of 2. On the other hand, Walch & Naab (2015) found that multiple SNe might increase the momentum per SN, by at least 25 per cent, depending on the delay time between SNe. The dependence on delay time further complicates this discrepancy between authors since neither set of authors used realistic delay time distributions for the number of SNe considered. Yadav et al. (2016) used 3D simulations to study how clustered SNe can merge into a superbubble (using a realistic SN delay time distribution), but did not study the momentum produced and did not test the effect of gas metallicity. So we are left with a series of questions: for a realistic delay time distribution of clustered SNe, does the momentum per SN increase or decrease relative to single SN models, and by how much? Does the result depend on the density or metallicity of the environment in which the SNe explode? Does it vary systematically with the number of SNe which are clustered together?

In this paper, we seek to measure directly the impact which clustering has on the momentum budget of SNe. In order to sample a wide range of densities, metallicities, and cluster sizes, we create a suite of several hundred 1D, spherically symmetric simulations. Using a 1D geometry means we lose the ability to simulate non-spherically symmetric inhomogeneities but in doing so we gain the ability to probe a far wider parameter space than any previous studies of multiple SNe, and to achieve far higher spatial resolutions than previous works studying momentum feedback. As we will show, both are necessary for understanding how clustering impacts momentum feedback.

The remainder of this paper is as follows. In Section 2, we discuss the numerical methods used in this study. The numerical results of our simulations are presented in Section 3. In Section 4, we use these results to build a model which can predict the momentum injection per SN as a function of density, metallicity, and number of SNe, in a form suitable for inclusion in subgrid and analytic models. We discuss the significance of our results and model in Section 5, comparing to previous works. Finally, we summarize our conclusions in Section 6.

2 NUMERICAL METHODS

Our simulations make use of a custom-built 1D spherically symmetric moving-mesh code which solves the finite-volume equations of compressible hydrodynamics. Our code includes radiative cooling and injection of mass and energy by both SNe and pre-SN winds. The star cluster is assumed to lie at the centre of our simulation, with our computational domain beginning just outside the cluster and extending outwards. SN ejecta and pre-SN winds are added to the innermost zone of this domain. We run these simulations until all SNe have occurred and the momentum reaches a maximum.

In the rest of this section, we go into greater depth on the numerical methods used in our simulations and the limitations of our assumptions. In Appendix A, we test our code against both the analytic Sedov solution and the earlier numerical results of Thornton et al. (1998) for isolated SNe, and verify that it reproduces them well. For the interested reader, our code has been publicly released.¹

2.1 Initial conditions

All our simulations begin with a star cluster of mass M_cluster placed at the origin surrounded by an initially uniform, stationary ideal gas with adiabatic index γ = 5/3. We vary M_cluster from 10²–10⁵ M_⊙ (using steps of 1 dex, with an additional step at 10^2.5 M_⊙ to better resolve a key region of parameter space). We explore gas mass densities in steps of 1 dex, ranging from ρ₀ = 1.33 × 10⁻³ to 10²m_H cm⁻³, where m_H = 1.67 × 10⁻²⁴ g is the mass of the hydrogen atom, corresponding to gas number densities of n₀ = ρ₀/(1.33m_H) H nuclei cm⁻³ for a helium mass fraction Y = 0.23 and a metal mass fraction Z = 0.02. We consider gas metallicities in steps of half a dex, ranging from Z₀ = 10⁻³ to 10^0.5 Z_⊙, excluding 10^−2.5 Z_⊙, where we have taken solar metallicity to be Z_⊙ = 0.02. The density and metallicity grids are chosen to closely match those explored by Thornton et al. (1998). Consistent with Thornton et al., we compute mean molecular weights assuming a fixed helium mass fraction, Y = 0.23, taking the remainder to be hydrogen (X = 1 − Z − Y). The gas has an initial temperature of 10⁴ K, but is allowed to cool to lower temperatures via radiation – see Section 2.3.

Our simulations start with 1024 zones, linearly spaced from radii R_in to R_out, which follow the scaling

\begin{eqnarray} R_\mathrm{out} = 300 \left( \frac{\rho _0}{1.33 m_\mathrm{H} {\rm cm}^{-3} } \right)^{-1/3} \left( \frac{M_\mathrm{cluster}}{100 \,\mathrm{M}_{\odot } } \right)^{1/3} \mathrm{pc} \end{eqnarray}

(1)

\begin{eqnarray} R_\mathrm{in} = 10^{-4} R_\mathrm{out}. \end{eqnarray}

(2)

This scaling is somewhat arbitrary and was set by initial tests; it was chosen to approximately reflect the final size of each simulation. If the outer boundary is too small, the domain will automatically extend when a shock nears the outer boundary.

2.2 Hydrodynamics

Our code solves the equations of compressible hydrodynamics in spherical symmetry using a moving-mesh finite-volume method, including source terms for radiative cooling. Our method is an extension of the one implemented by Duffell (2016). The equations we solve are

\begin{equation} \frac{ \mathrm{d} }{ \mathrm{d}t} \int \boldsymbol {U} \mathrm{d}V - \int \boldsymbol {F} \mathrm{d}A = \boldsymbol {S}, \end{equation}

(3)

where |$\boldsymbol {U} \mathrm{d} V$| is the vector of conserved quantities

\begin{equation} \boldsymbol {U} ={\left(\begin{array}{cccc}\begin{array}{l}\rho \\ \rho u_r \\ \rho e \\ \rho Z \end{array} \end{array}\right)} \end{equation}

(4)

for a density ρ, a bulk fluid velocity u, a specific total energy e, and a local metallicity Z. The quantity |$\boldsymbol {F}$| is the conservative flux, given by

\begin{equation} \boldsymbol {F} ={\left(\begin{array}{cccc}\begin{array}{lll}(u_r - w_r) \rho \\ (u_r - w_r) \rho u_r & +\, P \\ (u_r - w_r) \rho e & +\, P u_r + H \\ (u_r - w_r) \rho Z \\ \end{array} \end{array}\right)}, \end{equation}

(5)

where w is the computational mesh velocity, which is set to be the average velocity of the two zones adjacent to the boundary

\begin{equation} w_r^{(i+1/2)} = \frac{u_r^{(i)} + u_r^{(i+1)}}{2} \end{equation}

(6)

approximating Lagrangian hydrodynamics. Here P is the pressure given by

\begin{equation} P = (\gamma - 1) \rho e_\mathrm{int} \end{equation}

(7)

and e_int is the specific internal energy:

\begin{equation} e_\mathrm{int} = e - \frac{1}{2} v_r^2. \end{equation}

(8)

At the inner boundary, we enforce a zero flux boundary condition; the outer boundary condition does not matter, as we automatically add zones before the shock reaches the outer boundary (and we have assumed that the background is homogeneous).

This formulation implicitly introduces artificial mesh viscosity, particular at the inner boundary (i.e. wall heating), which leads to unphysically high temperatures. We counteract this by explicitly including an artificial conduction term, H. We use the artificial conduction prescription of Noh (1987, equation 2.3):

\begin{equation} H = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}h_0 \rho |\Delta u| \Delta e_\mathrm{int} + h_1 \rho c_\mathrm{s} \Delta e_\mathrm{int} & \Delta u < 0 \\ 0 & \Delta u > 0, \end{array}\right. \end{equation}

(9)

where c_s is the adiabatic sound speed, h₀ and h₁ are tunable constants, typically of order unity, and Δ represents the differential of a variable across adjacent zones. We chose these constants to be h₀ = 0 and h₁ = 0.1, which experimentation showed were the smallest values which were still sufficient to remove most unphysical wall heating. This parametrization is similar to physical conduction in the strong shock regime with a saturated conduction coefficient which has been lowered by a factor of a few by turbulence and magnetic fields (Cowie & McKee 1977).

In addition to these conservative fluxes, we also include non-conservative source terms

\begin{equation} \boldsymbol {S} = \boldsymbol {S}_\mathrm{hydro} + \boldsymbol {S}_\mathrm{cooling} + \boldsymbol {S}_\mathrm{winds} \end{equation}

(10)

\begin{equation} \boldsymbol {S}_\mathrm{hydro} ={\left(\begin{array}{cccc}0 \\ \int 2 (P/r)\, {\rm d}V \\ 0 \\ 0 \end{array}\right)} \end{equation}

(11)

\begin{equation} \boldsymbol {S}_\mathrm{cooling} ={\left(\begin{array}{cccc}0 \\ 0 \\ \dot{E}_\mathrm{cooling} \\ 0 \end{array}\right)} \end{equation}

(12)

\begin{equation} \boldsymbol {S}_\mathrm{winds} = \dot{M}_\mathrm{winds} \Delta t{\left(\begin{array}{ccc}\begin{array}{ll}1\\ u_\mathrm{winds} \\ (1/2)u^2_\mathrm{winds} + e_{\rm int, winds}\\ Z_\mathrm{winds} \end{array} \end{array}\right)}. \end{equation}

(13)

We defer a discussion of the cooling rate |$\dot{E}_\mathrm{cooling}$| to Section 2.3, and a discussion of the wind source term (which is only added to the innermost zone) to Section 2.5.1.

The conservative fluxes (excepting conduction) were solved using an HLLC Riemann solver (Toro, Spruce & Speares 1994) taken from the implementation of Duffell (2016). Artificial conduction and the non-conservative source terms were handled by operator splitting, solving each term individually.

By using a moving mesh with w_r ≈ u_r, we approximate Lagrangian hydrodynamics. This reduced numerical errors in the advective flux terms, as well as automatically adjusting to give higher resolution at locations with higher densities (assuming we start with a grid of uniform density). This improves our accuracy at shocks, where high densities lead to rapid cooling, which drives the subsequent evolution of the SNR. By using an approximately Lagrangian scheme, we can better resolve the dynamically important regions, without wasting computational time on the less important diffuse bubble.

For strong shocks, we need to set a limit on how much zones can expand or compress. The innermost zone (where SN energy is injected) will significantly expand, so we need to split it in order to retain accuracy where our blasts are being injected. For computational efficiency, we also need to allow zones to merge, because otherwise the zones near the shock become so thin that the computational cost of evolving them is prohibitive. We handle zone splitting and merging using the adaptive mesh algorithm implemented by Duffell (2016): zones thicker (thinner) than 10 (0.1) times the average zone thickness are split (merged).

To improve numerical stability in regions of highly supersonic flow, we also implement a dual-energy formalism. This approach counteracts the common problem in conservative, total energy codes such as ours that, at Mach numbers greater than unity, the internal energy e_int can be much smaller than the total energy e, so that small truncation errors in e can correspond to an order unity or larger error in e_int, and thus in the temperature and radiative cooling rate. Our dual-energy approach is as follows. For most zones and time-steps, we follow the updated procedure described above and derive e_int from the mass, momentum, and total energy (equation 8). However, in any zone and time-step where this procedure yields a value of e_int < 0, we instead compute the internal energy via

\begin{equation} e_\mathrm{int}(t + \Delta t) = e_\mathrm{int}(t) \left( \frac{{\rm d}V(t+\Delta t)}{{\rm d}V(t)}\right)^{1-\gamma } + \Delta e_\mathrm{cool}. \end{equation}

(14)

This includes adiabatic heating/cooling and radiative cooling; this ignores advective fluxes, which should be minimal for a pseudo-Lagrangian code, and conductive fluxes. The errors introduced by this dual-energy formalism have negligible effects on the overall dynamics and numerical conservation of energy.

2.3 Cooling

Cooling plays a significant role in SNR evolution, with most of the energy from the SN being radiated from a thin, dense shell. To include this cooling, we use the grackle chemistry and cooling library (Bryan et al. 2014; Kim et al. 2014), using operator splitting to evolve the thermal energy over each time-step. grackle sub-steps the thermal evolution using cooling rates pre-computed using cloudy (Ferland et al. 1998), assuming ionization equilibrium but not thermal equilibrium between metallicity-dependent optically thin cooling and a cosmological UV background at redshift z = 0 providing photoheating and photoionization (Haardt & Madau 2012). For simplicity, we only include heating from a cosmological background, rather than including galactic heating sources. We leave testing more realistic heating backgrounds and non-ionization equilibrium cooling models for a later work.

2.4 Cluster model

In order to test the SN momentum produced by a cluster, we need to determine the number of SNe from a cluster and when those SNe will occur. For each simulation with a given cluster mass, we use the slug2 code (da Silva, Fumagalli & Krumholz 2012, 2014; Krumholz et al. 2015) to draw the desired mass in stars from a Kroupa (2002) initial mass function (IMF), using the default ‘Stop-nearest’ policy. All stars above an initial mass of 8 M_⊙ are assumed to result in core-collapse SNe, after stellar lifetimes determined by the Geneva stellar evolution tracks assuming solar metallicity (Z = 0.014; Ekström et al. 2012). Generally, we find 1 SN per roughly 100 M_⊙ of stars, and those SNe occur roughly 3–40 Myr after the birth of the cluster. Given the power-law tail of the IMF, we expect most SN to come from relatively low mass stars, M_⋆ ≈ 8 M_⊙. We do not include Type Ia SNe for most of our simulations, but we do test the impact of short-delay Type Ia SNe in Section 5.3.2.

Since we are stochastically drawing an IMF, our results for low-mass clusters can depend significantly on the random seed. To minimize the uncertainty in our results induced by this stochasticity, we ran multiple realizations of the lowest cluster masses. Specifically, we ran nine realizations of each 10² M_⊙ cluster and four realizations of each 10^2.5 M_⊙ cluster.

This cluster model is not perfect. The stellar evolution tracks assume a single stellar metallicity, while ideally we would like the tracks to depend on the background metallicity for each simulation. Our cluster model also ignores the effects of stellar rotation and binarity. All of these can affect stellar lifetimes.

2.5 SN injection model

When an SN occurs, we add energy, mass, and metals to the innermost computational zone. For the energy, we adopt a fixed injection of 10⁵¹ erg per SN. For the mass and metallicity, we use the data of Woosley & Heger (2007), who provide a grid of SN mass and metal yields as a function of initial stellar mass over a range of initial masses from 12 to 120 M_⊙. Within this range, we linearly interpolate as a function of initial mass; outside this range, we use the nearest neighbour (i.e. stars with masses 8–12 M_⊙ are assumed to produce the same yield as 12 M_⊙ stars).

As with our cluster model, this SN model is imperfect. First, this model overpredicts the ejecta mass for stars with initial masses of 8–12 M_⊙ (the most common progenitors). For example, this model predicts that a 9 M_⊙ star will eject 9.4 M_⊙ of material. This is clearly unphysical, but the true ejected mass is comparable; Sukhbold et al. (2016) show that the true ejected mass (≈7.4 M_⊙) differs by less than 50 per cent from our simplified model. Overall, this will tend to overpredict the ejecta mass, biasing our results towards slightly more efficient cooling and slightly lower momenta.

Biasing our results in the opposite direction (for fixed cluster mass), we assume that all of our massive stars explode, even though some low-mass progenitors (8–9 M_⊙) may not explode (Woosley & Heger 2015) and some high-mass progenitors will collapse directly into black holes (Ertl et al. 2016; Sukhbold et al. 2016), a pathway which can depend significantly on the stellar metallicity (Pejcha & Thompson 2015). However, while models differ as to whether more massive stars explode, almost all models agree that all 9–12 M_⊙ stars should explode, so the total number of SNe should not change drastically.

More significantly, SN energies can vary with initial progenitor mass. In particular, Sukhbold et al. (2016) find that stars with initial masses of 9–12 M_⊙ explode with <0.7 × 10⁵¹ erg of energy. As these are the most common progenitors, this leads to an IMF-averaged explosion energy of ≈0.6–0.8 × 10⁵¹ erg, depending on the explosion model. While it is common to assume an explosion energy of 10⁵¹ erg (e.g. Thornton et al. 1998; Kim et al. 2014; Iffrig & Hennebelle 2015; Kim & Ostriker 2015; Martizzi et al. 2015; Walch & Naab 2015), in doing so we overpredict the average SN energy by a factor of 1.2–1.7.

As with our stellar evolution tracks, the data of Woosley & Heger (2007) are only for stars of solar metallicity, so we are unable to vary our SN model with the background metallicity. Moreover, we are combining the ejecta computed by Woosley & Heger (2007) with the lifetimes computed by Ekström et al. (2012); these make different assumptions about stellar evolution, and are not fully consistent. Theoretical uncertainties in stellar lifetimes are not too worrisome though; for low-mass clusters, we are dominated by stochastic scatter in the IMF; for high-mass clusters, the SNe are effectively a continuous wind. The models of Woosley & Heger (2007) and Ekström et al. (2012) also differ in the details of pre-SN mass-loss, but these discrepancies primarily exist for the most massive stars, which are the least common in our simulations.

2.5.1 Winds

If SN ejecta mass is important because more mass leads to faster cooling, then we cannot simply ignore pre-SN mass-loss; that mass has to go somewhere. The pre-SN mass-loss can be determined from the data of Woosley & Heger (2007), but that does not tell us when that mass was lost, or its physical properties when it was lost (e.g. metallicity, wind velocity, wind energy). For simplicity, we assume that pre-SN mass-loss occurs uniformly through a star's lifetime, as a wind with metallicity equal to the background metallicity, at a velocity of 10³ km s⁻¹, and a temperature of 10⁴ K. The total mass, metal mass, momentum, and energy of this wind are added to the innermost zone.

3 NUMERICAL RESULTS

We run a total of 672 simulations, sampling a three-dimensional parameter space composed of density ρ, metallicity Z, and cluster mass M_cluster. As an example of the outcome of our simulations, in Fig. 1 we show the density profile immediately after the last SN occurs in the simulation with ρ = 1.33 m_H cm⁻³, Z = Z_⊙, and M_cluster = 10⁵ M_⊙. In Fig. 2, we show the density profile for this simulation at the time when the radial momentum reaches its maximum, and we show the momentum as a function of time in Fig. 3 and the cumulative energy budget as a function of time in Fig. 4.

Figure 1.

Example density profile of a simulation with Z = Z_⊙, ρ = 1.33m_H cm⁻³, and M_cluster = 10⁵ M_⊙ (N_SNe = 1008), shortly after the last SN (t = 38 Myr).

Open in new tab Download slide

Figure 2.

The same simulation as the one shown in Fig. 1, except now at the moment of peak momentum (t = 285 Myr). The shock has weakened, causing it to thicken considerably.

Open in new tab Download slide

Figure 3.

The evolution of the momentum per SN of the cluster shown in Figs 1 and 2. The time of maximal momentum is marked by a vertical black dashed line; the duration of SN events is denoted by solid black ticks. For many SNe, the energy injection behaves more like a continuous wind rather than discrete explosions.

Open in new tab Download slide

Figure 4.

The evolution of the cumulative energy budget of the cluster shown in Figs 1 through 3. The time of maximal momentum is marked by a vertical black dashed line; SN times are denoted by solid black ticks. The kinetic component is measured directly from the simulation; the radiated component is inferred by comparing the decrease in total energy compared to the decrease which would have occurred if there were no SNe and the background cooled anyway; the thermal component is assumed to be whatever remains.

Open in new tab Download slide

In Table 1, we provide an overview of our results, extracting the following key parameters when all SNe have occurred and the momentum reaches a maximum: the peak momentum p; the time, t, at which the momentum reaches a maximum (defining t = 0 as the time of cluster formation); the radius of the shock, R, at this time (defined by the furthest zone with an overdensity compared to the background); the mass of the remnant, M_R, enclosed by the shock radius at this time; and the kinetic and internal energies, E_{R, kin} and E_{R, int}, enclosed by the shock radius at this time; finally, we also include a flag for untrustworthy results, which we explain in the next paragraph. In Table 2, we provide the time-dependent evolution of these parameters for every simulation and each snapshot before the time of peak momentum.

Table 1.

Open in new tab

Overview of numeric results. The table shown here is only a stub; it is provided in its entirety as a machine-readable table. The machine-readable table also includes a column with an id unique to each row, to allow cross-referencing with Table 2, which has been hidden here to save space.

ρ	Z	M_cluster	N_SNe	t	R	M_R	p	E_{R, kin}	E_{R, int}	Flag
(1.33 m_H cm⁻³)	(Z_⊙)	(M_⊙)		(Myr)	(pc)	( M_⊙)	(g cm s⁻¹)	(erg)	(erg)
10⁺⁰	10^+0.0	10^2.0	1	5.9	56.0	2.42 × 10⁴	5.22 × 10⁴³	3.14 × 10⁴⁹	1.77 × 10⁴⁸	0
10⁺⁰	10^+0.0	10^2.5	3	55.3	294.3	3.51 × 10⁶	1.10 × 10⁴⁵	9.47 × 10⁴⁹	2.73 × 10⁵⁰	0
10⁺⁰	10^+0.0	10^3.0	11	91.8	533.0	2.08 × 10⁷	6.72 × 10⁴⁵	5.94 × 10⁵⁰	1.63 × 10⁵¹	0
10⁺⁰	10^+0.0	10^4.0	104	173.3	1149.0	2.09 × 10⁸	7.18 × 10⁴⁶	6.67 × 10⁵¹	1.63 × 10⁵²	0
10⁺⁰	10^+0.0	10^5.0	1008	285.2	2133.0	1.34 × 10⁹	4.88 × 10⁴⁷	4.79 × 10⁵²	1.04 × 10⁵³	0

ρ	Z	M_cluster	N_SNe	t	R	M_R	p	E_{R, kin}	E_{R, int}	Flag
(1.33 m_H cm⁻³)	(Z_⊙)	(M_⊙)		(Myr)	(pc)	( M_⊙)	(g cm s⁻¹)	(erg)	(erg)
10⁺⁰	10^+0.0	10^2.0	1	5.9	56.0	2.42 × 10⁴	5.22 × 10⁴³	3.14 × 10⁴⁹	1.77 × 10⁴⁸	0
10⁺⁰	10^+0.0	10^2.5	3	55.3	294.3	3.51 × 10⁶	1.10 × 10⁴⁵	9.47 × 10⁴⁹	2.73 × 10⁵⁰	0
10⁺⁰	10^+0.0	10^3.0	11	91.8	533.0	2.08 × 10⁷	6.72 × 10⁴⁵	5.94 × 10⁵⁰	1.63 × 10⁵¹	0
10⁺⁰	10^+0.0	10^4.0	104	173.3	1149.0	2.09 × 10⁸	7.18 × 10⁴⁶	6.67 × 10⁵¹	1.63 × 10⁵²	0
10⁺⁰	10^+0.0	10^5.0	1008	285.2	2133.0	1.34 × 10⁹	4.88 × 10⁴⁷	4.79 × 10⁵²	1.04 × 10⁵³	0

Table 1.

Open in new tab

Overview of numeric results. The table shown here is only a stub; it is provided in its entirety as a machine-readable table. The machine-readable table also includes a column with an id unique to each row, to allow cross-referencing with Table 2, which has been hidden here to save space.

ρ	Z	M_cluster	N_SNe	t	R	M_R	p	E_{R, kin}	E_{R, int}	Flag
(1.33 m_H cm⁻³)	(Z_⊙)	(M_⊙)		(Myr)	(pc)	( M_⊙)	(g cm s⁻¹)	(erg)	(erg)
10⁺⁰	10^+0.0	10^2.0	1	5.9	56.0	2.42 × 10⁴	5.22 × 10⁴³	3.14 × 10⁴⁹	1.77 × 10⁴⁸	0
10⁺⁰	10^+0.0	10^2.5	3	55.3	294.3	3.51 × 10⁶	1.10 × 10⁴⁵	9.47 × 10⁴⁹	2.73 × 10⁵⁰	0
10⁺⁰	10^+0.0	10^3.0	11	91.8	533.0	2.08 × 10⁷	6.72 × 10⁴⁵	5.94 × 10⁵⁰	1.63 × 10⁵¹	0
10⁺⁰	10^+0.0	10^4.0	104	173.3	1149.0	2.09 × 10⁸	7.18 × 10⁴⁶	6.67 × 10⁵¹	1.63 × 10⁵²	0
10⁺⁰	10^+0.0	10^5.0	1008	285.2	2133.0	1.34 × 10⁹	4.88 × 10⁴⁷	4.79 × 10⁵²	1.04 × 10⁵³	0

ρ	Z	M_cluster	N_SNe	t	R	M_R	p	E_{R, kin}	E_{R, int}	Flag
(1.33 m_H cm⁻³)	(Z_⊙)	(M_⊙)		(Myr)	(pc)	( M_⊙)	(g cm s⁻¹)	(erg)	(erg)
10⁺⁰	10^+0.0	10^2.0	1	5.9	56.0	2.42 × 10⁴	5.22 × 10⁴³	3.14 × 10⁴⁹	1.77 × 10⁴⁸	0
10⁺⁰	10^+0.0	10^2.5	3	55.3	294.3	3.51 × 10⁶	1.10 × 10⁴⁵	9.47 × 10⁴⁹	2.73 × 10⁵⁰	0
10⁺⁰	10^+0.0	10^3.0	11	91.8	533.0	2.08 × 10⁷	6.72 × 10⁴⁵	5.94 × 10⁵⁰	1.63 × 10⁵¹	0
10⁺⁰	10^+0.0	10^4.0	104	173.3	1149.0	2.09 × 10⁸	7.18 × 10⁴⁶	6.67 × 10⁵¹	1.63 × 10⁵²	0
10⁺⁰	10^+0.0	10^5.0	1008	285.2	2133.0	1.34 × 10⁹	4.88 × 10⁴⁷	4.79 × 10⁵²	1.04 × 10⁵³	0

Table 2.

Open in new tab

Momentum evolution. The table shown here is only a stub; it is provided in its entirety as a machine-readable table.

ID	t	R_shock(t)	M_R(t)	p(t)	E_{R, kin}(t)	E_{R, int}(t)	N_SNe( < t)
	(Myr)	(pc)	( M_⊙)	(g cm s⁻¹)	(erg)	(erg)
25451948-485f-46fe-b87b-f4329d03b203	4.0	1.3	6.92 × 10⁰	0.00 × 10⁰	0.00 × 10⁰	1.00 × 10⁵¹	1
25451948-485f-46fe-b87b-f4329d03b203	5.0	61.4	3.20 × 10⁴	1.73 × 10⁴⁴	2.58 × 10⁵⁰	5.84 × 10⁵⁰	2
25451948-485f-46fe-b87b-f4329d03b203	10.2	152.9	4.92 × 10⁵	1.24 × 10⁴⁵	8.26 × 10⁵⁰	1.14 × 10⁵¹	5
25451948-485f-46fe-b87b-f4329d03b203	15.4	209.2	1.26 × 10⁶	2.18 × 10⁴⁵	9.88 × 10⁵⁰	1.40 × 10⁵¹	7
25451948-485f-46fe-b87b-f4329d03b203	20.5	249.3	2.13 × 10⁶	2.83 × 10⁴⁵	1.07 × 10⁵¹	2.43 × 10⁵¹	9
25451948-485f-46fe-b87b-f4329d03b203	25.6	283.1	3.12 × 10⁶	3.56 × 10⁴⁵	1.06 × 10⁵¹	1.61 × 10⁵¹	9
25451948-485f-46fe-b87b-f4329d03b203	35.3	335.1	5.18 × 10⁶	4.59 × 10⁴⁵	1.05 × 10⁵¹	2.84 × 10⁵¹	11

ID	t	R_shock(t)	M_R(t)	p(t)	E_{R, kin}(t)	E_{R, int}(t)	N_SNe( < t)
	(Myr)	(pc)	( M_⊙)	(g cm s⁻¹)	(erg)	(erg)
25451948-485f-46fe-b87b-f4329d03b203	4.0	1.3	6.92 × 10⁰	0.00 × 10⁰	0.00 × 10⁰	1.00 × 10⁵¹	1
25451948-485f-46fe-b87b-f4329d03b203	5.0	61.4	3.20 × 10⁴	1.73 × 10⁴⁴	2.58 × 10⁵⁰	5.84 × 10⁵⁰	2
25451948-485f-46fe-b87b-f4329d03b203	10.2	152.9	4.92 × 10⁵	1.24 × 10⁴⁵	8.26 × 10⁵⁰	1.14 × 10⁵¹	5
25451948-485f-46fe-b87b-f4329d03b203	15.4	209.2	1.26 × 10⁶	2.18 × 10⁴⁵	9.88 × 10⁵⁰	1.40 × 10⁵¹	7
25451948-485f-46fe-b87b-f4329d03b203	20.5	249.3	2.13 × 10⁶	2.83 × 10⁴⁵	1.07 × 10⁵¹	2.43 × 10⁵¹	9
25451948-485f-46fe-b87b-f4329d03b203	25.6	283.1	3.12 × 10⁶	3.56 × 10⁴⁵	1.06 × 10⁵¹	1.61 × 10⁵¹	9
25451948-485f-46fe-b87b-f4329d03b203	35.3	335.1	5.18 × 10⁶	4.59 × 10⁴⁵	1.05 × 10⁵¹	2.84 × 10⁵¹	11

Table 2.

Open in new tab

Momentum evolution. The table shown here is only a stub; it is provided in its entirety as a machine-readable table.

ID	t	R_shock(t)	M_R(t)	p(t)	E_{R, kin}(t)	E_{R, int}(t)	N_SNe( < t)
	(Myr)	(pc)	( M_⊙)	(g cm s⁻¹)	(erg)	(erg)
25451948-485f-46fe-b87b-f4329d03b203	4.0	1.3	6.92 × 10⁰	0.00 × 10⁰	0.00 × 10⁰	1.00 × 10⁵¹	1
25451948-485f-46fe-b87b-f4329d03b203	5.0	61.4	3.20 × 10⁴	1.73 × 10⁴⁴	2.58 × 10⁵⁰	5.84 × 10⁵⁰	2
25451948-485f-46fe-b87b-f4329d03b203	10.2	152.9	4.92 × 10⁵	1.24 × 10⁴⁵	8.26 × 10⁵⁰	1.14 × 10⁵¹	5
25451948-485f-46fe-b87b-f4329d03b203	15.4	209.2	1.26 × 10⁶	2.18 × 10⁴⁵	9.88 × 10⁵⁰	1.40 × 10⁵¹	7
25451948-485f-46fe-b87b-f4329d03b203	20.5	249.3	2.13 × 10⁶	2.83 × 10⁴⁵	1.07 × 10⁵¹	2.43 × 10⁵¹	9
25451948-485f-46fe-b87b-f4329d03b203	25.6	283.1	3.12 × 10⁶	3.56 × 10⁴⁵	1.06 × 10⁵¹	1.61 × 10⁵¹	9
25451948-485f-46fe-b87b-f4329d03b203	35.3	335.1	5.18 × 10⁶	4.59 × 10⁴⁵	1.05 × 10⁵¹	2.84 × 10⁵¹	11

ID	t	R_shock(t)	M_R(t)	p(t)	E_{R, kin}(t)	E_{R, int}(t)	N_SNe( < t)
	(Myr)	(pc)	( M_⊙)	(g cm s⁻¹)	(erg)	(erg)
25451948-485f-46fe-b87b-f4329d03b203	4.0	1.3	6.92 × 10⁰	0.00 × 10⁰	0.00 × 10⁰	1.00 × 10⁵¹	1
25451948-485f-46fe-b87b-f4329d03b203	5.0	61.4	3.20 × 10⁴	1.73 × 10⁴⁴	2.58 × 10⁵⁰	5.84 × 10⁵⁰	2
25451948-485f-46fe-b87b-f4329d03b203	10.2	152.9	4.92 × 10⁵	1.24 × 10⁴⁵	8.26 × 10⁵⁰	1.14 × 10⁵¹	5
25451948-485f-46fe-b87b-f4329d03b203	15.4	209.2	1.26 × 10⁶	2.18 × 10⁴⁵	9.88 × 10⁵⁰	1.40 × 10⁵¹	7
25451948-485f-46fe-b87b-f4329d03b203	20.5	249.3	2.13 × 10⁶	2.83 × 10⁴⁵	1.07 × 10⁵¹	2.43 × 10⁵¹	9
25451948-485f-46fe-b87b-f4329d03b203	25.6	283.1	3.12 × 10⁶	3.56 × 10⁴⁵	1.06 × 10⁵¹	1.61 × 10⁵¹	9
25451948-485f-46fe-b87b-f4329d03b203	35.3	335.1	5.18 × 10⁶	4.59 × 10⁴⁵	1.05 × 10⁵¹	2.84 × 10⁵¹	11

Not all of our simulations are trustworthy. In some runs, a strong reverse shock reaches the inner boundary before the SNR momentum peaks. In these simulations, the shock reflects off our hard inner boundary, whereas in reality the shock converging on the origin would certainly become unstable and would not reflect. In these cases, we cannot reasonably measure a maximum momentum. Fortunately, this behaviour only occurs in a small part of our parameter space (40/672 runs), and in what follows we will exclude these runs from our analysis. We also exclude any other realizations of the same initial conditions (an additional 99/672 runs) so as not to bias ourselves by only allowing atypical realizations. In Section 4.1.1, we explain the astrophysical causes and implications of these flagged runs.

The quantities M_R and E_{R, int} need to be interpreted with some care. At late times (when these quantities are extracted), the shock has weakened and is becoming a linear sound wave moving through a uniform medium, as illustrated in Fig. 2. At this time, the bubble-shell decomposition no longer is a good description, and M_R and E_{R, int} are becoming increasingly dominated by background material which has simply had a sound wave pass through it, but has not been irreversibly affected by a shock. It would be inaccurate to include this material and its internal energy in SN-driven ‘feedback’, and it is difficult to meaningfully disentangle SN-dominated material and background-dominated material as the SNR is merging into the ISM. It is easier to disentangle kinetic variables since the background is static; all momentum and kinetic energy must be a result of the SNe. The kinetic energy E_{R, kin} does not asymptote, but it varies slowly at late times, as illustrated in Fig. 4.

As an example of how these results vary across our parameter space, we plot the asymptotic momentum per SN in two cuts through this parameter space in Figs 5 (momentum per SN as a function of ρ and N_SNe at fixed Z) and 6 (momentum per SN as a function of Z and N_SNe at fixed ρ). In Figs 7 and 8, we provide analogous figures for the asymptotic kinetic energy, which is typically 1–10 per cent of the injected SN energy.

Figure 5.

An overview of the final momentum per SN and how it varies with the number of SNe and gas density at fixed metallicity (Z = Z_⊙). The locations of our simulations in parameter space are marked by black scatter points (excluding flagged runs), which are not a perfect grid because the numbers of SNe are drawn stochastically. The colour image is an interpolation of our simulation results using a Gaussian radial basis function, evaluated on a 100 by 100 grid with seven grey-scale contours logarithmically spaced between 3 × 10³ and 6 × 10⁴ km s⁻¹ (exclusive).

Open in new tab Download slide

The same as Fig. 5, except now allowing metallicity to vary while holding density fixed at 1.33 mH cm−3. The top contour level shown in Fig. 5 is not shown here, as the dynamic range of the data is not as large.

Figure 6.

The same as Fig. 5, except now allowing metallicity to vary while holding density fixed at 1.33 m_H cm⁻³. The top contour level shown in Fig. 5 is not shown here, as the dynamic range of the data is not as large.

Open in new tab Download slide

Same as Fig. 5, except now the colour image shows final kinetic energy with four grey-scale contours linearly spaced between 2 × 1049 and 1.2 × 1050 erg (exclusive).

Figure 7.

Same as Fig. 5, except now the colour image shows final kinetic energy with four grey-scale contours linearly spaced between 2 × 10⁴⁹ and 1.2 × 10⁵⁰ erg (exclusive).

Open in new tab Download slide

The same as Fig. 7, except now allowing metallicity to vary while holding density fixed at 1.33 mH cm−3. The top two contour levels shown in Fig. 7 are not shown here, as the dynamic range of the data is not as large.

Figure 8.

The same as Fig. 7, except now allowing metallicity to vary while holding density fixed at 1.33 m_H cm⁻³. The top two contour levels shown in Fig. 7 are not shown here, as the dynamic range of the data is not as large.

Open in new tab Download slide

4 THE MOMENTUM BUDGET OF CLUSTERED SNe

Figs 5 and 6 show significant structure in the momentum as a function of number of SNe, gas metallicity, and density. In particular, we note three behaviours. (1) For fixed density and metallicity, starting with a few SNe, the momentum per SN initially increases with increasing number of SNe, reaches a maximum between 10–100 SNe (the exact location depends on density and metallicity), then decreases. (2) For a few SNe, the momentum per SN increases with decreasing density and gas metallicity. (3) For many SNe, the opposite is true, as momentum per SN increases with increasing density and, to a smaller extent, metallicity. In this section, we show how these primary behaviours are a consequence of clustered SNR evolution falling into one of two physical regimes: the small-N regime and the superbubble regime.

4.1 Qualitative analysis

4.1.1 The small-N regime

To understand how SN momentum budgets act when the number of clustered SNe is relatively small, we start with its limiting case: single, isolated SNe. Feedback from isolated SNe has been well explored, as discussed in Section 1. In particular, Thornton et al. (1998) found that lower gas metallicities and densities resulted in higher energy feedback, which is what we expect physically; lower gas metallicity and density result in weaker cooling, sapping less energy from an SNR, increasing the amount of energy feedback. This is also expected to apply for momentum feedback, and in Fig. 9 we show that our results match the scaling between momentum and density expected by Cioffi et al. (1988).

Figure 9.

The scaling of momentum with background density, for Z = Z_⊙ and N_SNe = 1, compared to the p ∝ ρ^−1/7 scaling (normalized to the mean momentum of the lowest density clusters) expected for isolated SNe in a homogeneous background (Cioffi et al. 1988).

Open in new tab Download slide

As the number of SNe increases, the picture is similar to a series of isolated SNe, but with each successive SN occurring in a lower density bubble. As discussed above, this leads to progressively more efficient momentum production as the region is progressively evacuated. Fig. 10 illustrates this process directly, by plotting the momentum versus time for a simulation in which two SNe occur. The first SN occurs 20 Myr after cluster formation and its remnant quickly asymptotes to a momentum ≈3 × 10⁵ M_⊙ km s⁻¹, in agreement with the usual value found for single SNe. The second SN occurs 5 Myr later, and thanks to the vastly lower density inside the bubble, experiences much smaller radiative losses. This leads it to inject ≈2 × 10⁶ M_⊙ km s⁻¹ of momentum by the time the momentum peaks, which is almost 10 times more momentum than injected by the first SN.

The evolution of the momentum per SN of a Z = Z⊙, ρ = 1.33 mH cm−3, and NSNe = 2 cluster. The moment of maximal momentum is marked by the vertical dashed black line; the times of SNe are denoted by solid black ticks.

Figure 10.

The evolution of the momentum per SN of a Z = Z_⊙, ρ = 1.33 m_H cm⁻³, and N_SNe = 2 cluster. The moment of maximal momentum is marked by the vertical dashed black line; the times of SNe are denoted by solid black ticks.

Open in new tab Download slide

This breaks down for the clusters with the fewest SNe embedded in the highest density backgrounds, which behave more like multiple, isolated SNe. As density increases, SNRs evolve more rapidly, quickly cooling, lowering their internal pressure, and then being crushed by the pressure of the surrounding ISM. For the most dense gas with the fewest SNe, the SNR bubbles can be destroyed between SNe. As subsequent SNe are no longer occurring in the bubble of previous SNe, we believe that clustering effects should be minimal.

Unfortunately, since the bubble has collapsed before all SNe have been injected, our numerical methods break down due to the reverse shock propagating all the way to the origin and undergoing unphysical reflection (see Section 3). We are therefore forced to exclude this regime from our analysis. In our three-dimensional parameter space, the excluded region is roughly defined by the parameters

\begin{eqnarray*} \rho \ge 1.33 \, m_\mathrm{H}\, {\rm cm}^{-3}\quad \mathrm{and}\quad N_\mathrm{SNe} < 10, \quad\mathrm{and}\quad Z < 0.1 \,\mathrm{Z}_{\odot }. \end{eqnarray*}

While we cannot simulate this part of parameter space directly, the above analysis suggests that there should be no clustering effects present in it, and thus for the purpose of subgrid modelling it is likely safe to adopt a momentum budget of 3 × 10⁵ M_⊙ km s⁻¹ per SN, the same as in the isolated SN regime.

4.1.2 The superbubble regime

While the few-SN model predicts that momentum efficiency increases as the number of SNe increases, our data show a turnover after about 10–100 SNe, beyond which the momentum efficiency begins to drop as the number of SNe increases. This can be understood within the framework of a superbubble powered by a continuous wind (see Fig. 3 for an example of the momentum evolution of a large cluster). As more SNe occur, the bubble density decreases while the bubble temperature increases, both of which lead to less efficient cooling. While this leads to strong momentum feedback for a few SNe, it eventually saturates; if most of the energy is already being retained, suppressing cooling even further will only have a marginal effect.

Castor, McCray & Weaver (1975) provide a simplified bubble model which allows us to begin to understand superbubble evolution. They assume a constant energy injection rate, but if that energy is injected over a period of time which is the same for all clusters (effectively assuming that stellar evolution models do not depend strongly on metallicity or density), then their approach is also valid for an energy injection rate which is a power law with respect to time. Using their bubble solution, we can find the momentum per SN at the time of the last SN:

\begin{equation} p(t_{\rm last SN})/N_\mathrm{SNe} \propto N_\mathrm{SNe}^{-0.2} \rho _0^{0.2}. \end{equation}

(15)

We compare this predicted scaling to our numerical results for the lowest density simulated (thereby ensuring that we are as close as possible to the adiabatic limit) in Fig. 11. As the plot shows, the analytic scaling is in reasonable agreement with the numeric results.

Figure 11.

The scaling of momentum per SN with number of SNe, evaluated at the time of the last SN for each cluster. The clusters shown all have solar metallicity and the lowest density simulated, 1.33 × 10⁻³m_H cm⁻³. We plot the theoretical scaling for an adiabatic superbubble (equation 15), normalized to the cluster with the most SNe (the cluster which is expected to best correspond to the adiabatic case).

Open in new tab Download slide

As shown in Fig. 3, a significant amount of momentum evolution occurs after the last SN. During this phase, the superbubble expands adiabatically until the bubble pressure equals the ISM pressure, P₀, at which point the shell's momentum reaches a maximum, since the pressure gradient switches direction. For an adiabatic index γ = 5/3 for the gas inside the SNR, this results in a final momentum per SN

\begin{eqnarray} p_\mathrm{final}/N_\mathrm{SNe} &\propto & P_0^{-1/(2 \gamma )} N_\mathrm{SNe}^{-0.2 + (0.2/\gamma )} \rho _0^{0.2 + (0.3/\gamma )} \nonumber \\ &\approx & P_0^{-0.3} N_\mathrm{SNe}^{-0.08} \rho _0^{0.38}. \end{eqnarray}

(16)

This analytic scaling with respect to number of SNe can be compared to our numeric data in Fig. 12; the scaling with respect to gas density is shown in Fig. 13.

The scaling of asymptotic momentum per SN with number of SNe. These are the same clusters as those shown in Fig. 11 (1.33 × 10−3mH cm−3; Z = Z⊙) but evaluated at a different time. We plot theoretical scalings for an adiabatic superbubble at time of the last SN (blue solid line; equation 15) and at the time the interior pressure equals the exterior pressure (green dashed line; equation 16).

Figure 12.

The scaling of asymptotic momentum per SN with number of SNe. These are the same clusters as those shown in Fig. 11 (1.33 × 10⁻³m_H cm⁻³; Z = Z_⊙) but evaluated at a different time. We plot theoretical scalings for an adiabatic superbubble at time of the last SN (blue solid line; equation 15) and at the time the interior pressure equals the exterior pressure (green dashed line; equation 16).

Open in new tab Download slide

The scaling of asymptotic momentum per SN with background density, for Z = Z⊙ and Mcluster = 105 M⊙ (NSNe ≈ 103) clusters, compared to the superbubble predictions for the end of the SN injection phase (p/NSNe ∝ ρ0.2) and when the interior pressure equals the exterior pressure (p/NSNe ∝ ρ0.2 + (0.3/γ)).

Figure 13.

The scaling of asymptotic momentum per SN with background density, for Z = Z_⊙ and M_cluster = 10⁵ M_⊙ (N_SNe ≈ 10³) clusters, compared to the superbubble predictions for the end of the SN injection phase (p/N_SNe ∝ ρ^0.2) and when the interior pressure equals the exterior pressure (p/N_SNe ∝ ρ^{0.2 + (0.3/γ)}).

Open in new tab Download slide

It is not surprising that these scalings are not perfect; Sharma et al. (2014) predict that even 10³ SNe are not enough to satisfy assumptions behind models like those of Castor et al. (1975). Specifically, Sharma et al. (2014) predict no wind-dominated region (where ρ ∝ r⁻²) which ends in a stable termination shock before the pressure-dominated bubble begins; these predictions are in agreement with our results (see Figs 1 and 2). Our simulations do not satisfy all of the assumptions of superbubble models like those of Castor et al. (1975); these models are sufficient for a qualitative analysis, but they are insufficient for a quantitative understanding.

4.2 Quantitative model

Informed by the qualitative understanding developed in Section 4.1, we now construct a quantitative parametric model which we constrain using our simulation results (Table 1). In both the small-N and superbubble limits, we expect the results to behave like a power law in number of SNe, gas density, and metallicity, but we expect these to be different power laws. Therefore, we choose to construct a model of two power laws with a smooth break. In the few-SN (small-N) limit, we use a model of the form

\begin{eqnarray} \left(\frac{p}{N_\mathrm{SNe}}\right)_\mathrm{few} &=& \left( \frac{p}{N_\mathrm{SNe}} \right)_{0,\mathrm{few}} \left(\frac{Z}{\mathrm{Z}_{\odot }}\right)^{\eta _{Z,\mathrm{few}}} \nonumber\\ &&\times \left(\frac{\rho }{m_\mathrm{H} \, \mathrm{cm}^{-3}}\right)^{\eta _{\rho ,\mathrm{few}}} \left(\frac{N_\mathrm{SNe}}{1}\right)^{\eta _{N,\mathrm{few}}} \end{eqnarray}

(17)

and in the many SNe (superbubble) limit, we use a similar form,

\begin{eqnarray} \left(\frac{p}{N_\mathrm{SNe}}\right)_\mathrm{many} &=& \left( \frac{p}{N_\mathrm{SNe}} \right)_{0,\mathrm{many}} \left(\frac{Z}{\mathrm{Z}_{\odot }}\right)^{\eta _{Z,\mathrm{many}}} \nonumber\\ &&\times \left(\frac{\rho }{m_\mathrm{H} \, \mathrm{cm}^{-3}}\right)^{\eta _{\rho ,\mathrm{many}}} \left(\frac{N_\mathrm{SNe}}{1000}\right)^{\eta _{N,\mathrm{many}}} \end{eqnarray}

(18)

which are smoothly combined using

\begin{eqnarray} \frac{p}{N_\mathrm{SNe}} = \frac{\left(\frac{p}{N_\mathrm{SNe}}\right)_\mathrm{few} \left(\frac{p}{N_\mathrm{SNe}}\right)_\mathrm{many}}{\left(\frac{p}{N_\mathrm{SNe}}\right)_\mathrm{few}+\left(\frac{p}{N_\mathrm{SNe}}\right)_\mathrm{many}} \end{eqnarray}

(19)

\begin{eqnarray} \qquad\,\approx \mathrm{min}\left[ \left(\frac{p}{N_\mathrm{SNe}}\right)_\mathrm{few},\left(\frac{p}{N_\mathrm{SNe}}\right)_\mathrm{many} \right]. \end{eqnarray}

(20)

Assuming that our simulation results have a random additive Gaussian noise of variance σ², we can construct a Gaussian likelihood function for the results of each simulation. Even though σ² is unknown, this allows us to determine a maximum likelihood estimate (MLE) for our best-fitting model parameters. We would also like to understand the uncertainties in those parameters. For that we need the posterior, which through Bayes’ theorem requires a prior distribution, π on those parameters. Since we do not have strong prior information on most of these parameters, we choose uniform, independent priors on our parameters: log (σ²), log (p/N_SNe)_{0, few}, η_{Z, few}, η_{ρ, few}, η_{N, few}, log (p/N_SNe)_{0, many}, η_{Z, many}, η_{ρ, many}, η_{N, many}.

Combining this prior with a Gaussian likelihood for our data results in the posterior distribution of our model parameters. We sample this posterior distribution using a Markov chain Monte Carlo scheme, using Gibbs sampling to draw samples of σ² from an inverse gamma distribution and using a Metropolis–Hastings random walk for the remaining parameters. We use the MLE as the starting guess, discard the first 10 000 steps as burn-in steps, and save the next 100 000 steps. Using these samples, we can now estimate uncertainties on our model parameters: for each parameter, we use the median as our best-fitting value, and the 16th and 84th percentiles as our uncertainty interval (effectively marginalizing over all other parameters) resulting in

\begin{eqnarray} \left( \frac{p}{N_\mathrm{SNe}} \right)_{0,\mathrm{few}} = 4249^{+741}_{-683} \times 100 \, \,\mathrm{M}_{\odot }\, {\rm km\, s}^{-1} \end{eqnarray}

(21)

\begin{eqnarray} \eta _{Z,\mathrm{few}} = 0.05^{+0.05}_{-0.06} \end{eqnarray}

(22)

\begin{eqnarray} \eta _{\rho ,\mathrm{few}} = -0.06^{+0.03}_{-0.03} \end{eqnarray}

(23)

\begin{eqnarray} \eta _{N,\mathrm{few}} = 2.20^{+0.24}_{-0.23} \end{eqnarray}

(24)

\begin{eqnarray} \left( \frac{p}{N_\mathrm{SNe}} \right)_{0,\mathrm{many}} &= 23546^{+1072}_{-1073} \times 100 \, \,\mathrm{M}_{\odot }\, {\rm km\, s}^{-1} \end{eqnarray}

(25)

\begin{eqnarray} \eta _{Z,\mathrm{many}} = 0.15^{+0.01}_{-0.01} \end{eqnarray}

(26)

\begin{eqnarray} \eta _{\rho ,\mathrm{many}} = 0.14^{+0.01}_{-0.01} \end{eqnarray}

(27)

\begin{eqnarray} \eta _{N,\mathrm{many}} = -0.07^{+0.02}_{-0.02} \end{eqnarray}

(28)

\begin{eqnarray} \sigma = 6075^{+214}_{-202} \times 100 \, \,\mathrm{M}_{\odot } \,{\rm km\, s}^{-1} \end{eqnarray}

(29)

Our posterior samples are also useful for estimating the uncertainty in the predicted momentum from a particular cluster. For a given gas metallicity, density, and number of SNe, each posterior sample predicts a slightly different momentum and an uncertainty σ on that momentum. For each posterior sample, we then generate realizations of the noise with variance σ². For N posterior samples and M noise realizations, this gives us N × M samples of the momentum predictive distribution, which allows us to estimate the median and the 16th and 84th percentiles of the predictive distribution for a given cluster. We compare this predictive model and its uncertainties to our a subset of our numeric data in Fig. 14.

Comparison of a slice of our simulation results (Z = Z⊙, ρ = 1.33mH cm−3) to our model with an uncertainty envelop bounding the 16th and 84th percentiles of our predictive momentum model. Some slices fit better and some slices fit worse, but overall this is a representative slice. For comparison, we also plot a typical unclustered model, p/(100 M⊙NSNe) = 3000 km s−1 (Ostriker & Shetty 2011, green dashed line).

Figure 14.

Comparison of a slice of our simulation results (Z = Z_⊙, ρ = 1.33m_H cm⁻³) to our model with an uncertainty envelop bounding the 16th and 84th percentiles of our predictive momentum model. Some slices fit better and some slices fit worse, but overall this is a representative slice. For comparison, we also plot a typical unclustered model, p/(100 M_⊙N_SNe) = 3000 km s⁻¹ (Ostriker & Shetty 2011, green dashed line).

Open in new tab Download slide

5 DISCUSSION

Using our numeric results and quantitative model, we can now comment on the significance of these results in the context of previous works. We first examine the implications of our results for models of momentum-regulated star and galaxy formation in Section 5.1. We then compare our results to those of previous authors in Section 5.2, and in Section 5.3 we discuss the potential importance of physical processes we have omitted.

5.1 Implications of high momentum efficiency of clustered SNe

In models of momentum-driven feedback, the key parameter is p/m_*, the amount of momentum injected per unit mass of stars formed in a given system. Non-clustered models of SN momentum production usually assume |$p/m_{\ast } \approx p/(100 \,\mathrm{M}_{\odot } N_\mathrm{SNe}) \approx 1000 {\rm -} 3000$| km s⁻¹ (with a weak dependence on density) for a mass m_* of stars formed (Thompson, Quataert & Murray 2005; Ostriker & Shetty 2011; Shetty & Ostriker 2012; Dekel & Krumholz 2013; Faucher-Giguère, Quataert & Hopkins 2013; Hopkins et al. 2014; Hayward & Hopkins 2016; Kim & Ostriker 2015; Kimm et al. 2015). For a star cluster with a single SN (M_cluster ≈ 100 M_⊙), our best-fitting model is a little higher than but still consistent with 1000–3000 km s⁻¹ given the uncertainties in our model. For higher mass clusters, the discrepancy becomes significant. The most extreme difference is found for M_cluster = 10³–10⁴ M_⊙ (N_SNe = 10¹–10²), for which our value of p/m_* can be greater than the unclustered value by an order of magnitude (see Figs 5 and 14). But these are just the extremes; for a typical distribution of cluster masses found in a galaxy, what is the average effect?

To evaluate the mean value of p/m_* for star formation on galactic scales, we must integrate our model for individual clusters, p/M_cluster, over a cluster mass function dN/dM_cluster. The resulting mean momentum yield per unit mass of stars formed is

\begin{equation} \frac{p}{m_{\ast }} = \int \left(\frac{p}{M_\mathrm{cluster}}\right) \frac{{\rm d}N}{{\rm d} \ln M_{\rm cluster}} \, {\rm d}M_{\rm cluster}. \end{equation}

(30)

If we adopt a typical mass distribution |${\rm d}N/{\rm d}M_{\rm cluster} \propto M_{\rm cluster}^{-2}$| over the range M_cluster = 10²–10⁵ M_⊙, comparable to what is observed in nearby galaxies (Krumholz 2014, and references therein), and use our fitting formula (equation 19) to evaluate p/N_SNe as a function of M_cluster (for N_SNe ≈ M_cluster/100 M_⊙), this yields a value of p/m_* ≈ 1–2 × 10⁴ km s⁻¹ over the metallicity range Z/Z_⊙ = 0.01–1 and density range ρ/m_H = 0.1–10⁵. This is ∼0.5–1 dex higher than the value usually adopted based on single SN models. This result is only logarithmically sensitive to the adopted limits on the cluster mass function.

This increased momentum yield will significantly alter the conclusions of analytic models in which star formation is regulated primarily by SN momentum input (e.g. Ostriker & Shetty 2011; Shetty & Ostriker 2012; Faucher-Giguère et al. 2013; Hayward & Hopkins 2016). The same is true for models where SN momentum is primarily responsible for launching galactic winds (e.g. Dekel & Krumholz 2013; Hayward & Hopkins 2016; Thompson & Krumholz 2016). In general, the higher momentum yield we obtain will shift such models to predict lower star formation rates for fixed galactic surface densities, which may require re-tuning of other parameters to bring the models back into agreement with observed relationships between star formation rate and gas content.² The models will also predict stronger outflows, though these are significantly less constrained by observations.

The momentum yield per SN is also a critical input to numerical methods which handle subgrid feedback through explicit momentum injection (e.g. Kim et al. 2011; Hopkins et al. 2014; Kimm et al. 2015; Goldbaum et al. 2016). These models should be also be rerun using our updated estimates of the SN momentum yield. Even models which do not use explicit momentum injection, but which attempt to include SN feedback by explicitly resolving the Sedov phase (e.g. Hopkins et al. 2011), may need to be reconsidered, at least for simulations of galaxies large enough for there to be significant numbers of clustered SNe.

5.2 Comparison to previous work and convergence study

We find that clustered SNe generally lead to an increase in the momentum injected by SNe, in some cases by an order of magnitude. The results of previous authors diverge strongly, as noted in Section 1, with some finding that clustering leads to an enhancement in momentum per SN and others finding a decrease, but none finding an increase as large as an order of magnitude. To understand this discrepancy, we need to understand the role of mixing and how it enters various simulations.

In SN-driven bubbles, the cooling rate plays a significant role in setting the dynamics of the system. This cooling rate is itself affected by the mixing rate at interfaces between hot diffuse gas (which dominates the thermal energy) and cold dense gas (which is most radiatively efficient). If the mixing of energy and matter increases, the cooling rate can increase and the final momentum can decrease. This mixing can be increased by both non-physical sources (i.e. due to numerical diffusion) and physical sources (i.e. conduction or hydrodynamic instabilities which transport energy across the contact discontinuity). While these two channels have very different causes, they can have similar effects on the momentum and energy evolution of the system (Fierlinger et al. 2016). We will look at these two channels in turn.

Hydrodynamic solvers which advect mass between adjacent cells fundamentally introduce mixing errors. These errors can be decreased by improving the resolution or decreasing the mass advected between cells (as we have done by moving our numerical mesh with the fluid). In order to understand how the choice of resolution and numerical methods affects our results, we re-ran one of our clusters [ρ = 1.33 m_H cm⁻³, Z = Z_⊙, M_cluster = 10³ M_⊙, N_SNe = 11; this particular cluster will be useful in later comparisons to Kim & Ostriker (2015)] at a number of initial resolutions, using both a fixed mesh and a moving mesh. (A fixed mesh corresponds to Eulerian hydrodynamics – fixing w = 0 in our methods discussed in Section 2 – which is less accurate and more diffusive than our pseudo-Lagrangian methods.) We show the results in Fig. 15 and note a few key observations. First, the Lagrangian runs appear to be converged (within the uncertainties of our predictive model) by the resolution used for this cluster in our parameter study (an initial resolution of 0.6 pc). Secondly, the Eulerian runs introduce larger errors (as expected) and converge much more slowly. This suggests that Eulerian and low-resolution simulations could have greater errors than high-resolution, Lagrangian simulations. We can better understand these errors by comparing the energy evolution of a high-resolution Lagrangian run (Fig. 16) and a low-resolution Eulerian run (Fig. 17). While the same amount of energy is injected for each SN in both simulations, that energy is radiated away much more rapidly in the low-resolution (Eulerian) simulation, draining the bubble of the energy which drives the momentum growth seen in Fig. 15. [This connection between cooling time and resolution for multiple SNe was also found by Krause et al. (2013).] The amount of mixing can significantly impact the final momentum, but given the convergence seen in Fig. 15, the results we have obtained appear to be converged.

The momentum evolution of a ρ = 1.33 mH cm−3, Z = Z⊙, Mcluster = 103 M⊙ (NSNe = 11) cluster, rerun with a range of initial resolutions, using both Eulerian and Lagrangian methods. The asymptotic momentum predicted by our model is shown by the blue horizontal band which bounds the 16th and 84th percentiles of the predictive distribution.

Figure 15.

The momentum evolution of a ρ = 1.33 m_H cm⁻³, Z = Z_⊙, M_cluster = 10³ M_⊙ (N_SNe = 11) cluster, rerun with a range of initial resolutions, using both Eulerian and Lagrangian methods. The asymptotic momentum predicted by our model is shown by the blue horizontal band which bounds the 16th and 84th percentiles of the predictive distribution.

Open in new tab Download slide

Same as Fig. 4, except now for a ρ = 1.33 mH cm−3, Z = Z⊙, Mcluster = 103 M⊙ (NSNe = 11) cluster, evolved using Lagrangian methods with an initial resolution of 0.6 pc.

Figure 16.

Same as Fig. 4, except now for a ρ = 1.33 m_H cm⁻³, Z = Z_⊙, M_cluster = 10³ M_⊙ (N_SNe = 11) cluster, evolved using Lagrangian methods with an initial resolution of 0.6 pc.

Open in new tab Download slide

Figure 17.

Same as Fig. 16, except with the resolution degraded to 2.5 pc, and using a fixed Eulerian mesh.

Open in new tab Download slide

While our resolution study suggests that high-resolution, low-diffusion simulations are required to achieve accurate results, that conclusion might not apply if there are stronger, physical diffusive processes present (Fierlinger et al. 2016). For an SNR or a superbubble, a number of processes can mix gas between the hot bubble interior and the cool shell; for a review, see appendix B of Fierlinger et al. (2016). Our 1D code cannot simulate many of these processes directly,³ but higher dimensional simulations can. In order to test the effects of these more complex mixing interfaces, we will compare our results to existing 3D simulations of multiple SNe in an inhomogeneous background.

Martizzi et al. (2015), Walch & Naab (2015), and Kim & Ostriker (2015) all tested the effects which a turbulent or multiphase background might have on SNR evolution. For 1 SN, they all found that an inhomogeneous background makes a relatively small difference: a change of less than 60 per cent. They also test multiple SNe in an inhomogeneous background, but none compare the results to multiple SNe in a homogeneous background. So in order to understand the effect of mixing in the case of multiple SNe, we will compare one of our clusters (ρ = 1.33 m_H cm⁻³, Z = Z_⊙, M_cluster = 10³ M_⊙, N_SNe = 11; the cluster used in our resolution study) with a multiphase multiple SN cluster from Kim & Ostriker (2015, ρ = 1.4 m_H cm⁻³, Z = Z_⊙, N_SNe = 10). Note that although these clusters were chosen to be as similar as possible (except with a difference in background media), they also differ in SN delay time distributions, ejecta, and mass-loss prescriptions. Nevertheless, we can compare our results to those of Kim & Ostriker (2015) and find that their cluster cools much more rapidly than ours, leading to an asymptotic momentum per SN which is a factor of 20 lower than ours.

We believe that this difference is due to physical mixing which is present in their simulations but not ours. While our results can be brought into agreement with theirs by degrading our resolution and using less accurate numerical methods (as shown by our convergence study), subsequent simulations have shown their results to be converged with respect to resolution (Kim, Ostriker & Raileanu 2016). If a two-phase background results in significantly increased physical mixing, that would explain how their converged simulations could appear similar to our less accurate, unresolved simulations: a strong source of physical mixing can appear similar to strong artificial mixing (Fierlinger et al. 2016), while being less sensitive to the resolution.

One might ask at this point whether the enhanced momentum injection we find is solely a result of our use of 1D simulations, which necessarily suppress mixing. Such a conclusion might be comforting, but is far from warranted. The overall lesson to draw from this comparison is that the cooling rate and momentum budget for bubbles produced by multiple SNe is exquisitely sensitive to the amount of mixing, whether physical or numerical. For a homogeneous background, very high resolution is required to get a converged value for the asymptotic momentum. This resolution requirement can be rendered irrelevant if strong, physical mixing occurs, allowing convergence at much lower resolution. But accurate results require more than just convergence; to be confident in the accuracy of a set of results, one must be confident that the physical mixing processes have been properly captured. In their multiple SN simulations, Kim & Ostriker (2015) include a two-phase background, but not magnetic fields, which are known to suppress mixing across contact discontinuities in other contexts (e.g. Markevitch & Vikhlinin 2007). Thus, their results should probably be regarded as lower limits. Sharma et al. (2014) do include magnetic fields, but only the context of a uniform medium. Given the state of the field, and the resolution requirements we have obtained in the uniform case, it seems clear that there is an urgent need to re-examine the momentum budget of clustered SNe in a multidimensional context, properly including all the mechanisms which can both enhance and suppress mixing.

5.3 The effects of additional physics

In order to render the problem as clean as possible, we have focused only on Type II SN feedback in a uniform medium. We now consider how other physical processes which we have heretofore neglected might alter our results.

5.3.1 Pre-SN radiative feedback

Before any SNe occur, we expect pre-SN feedback to already be sculpting the region. In particular, ionizing radiation from young stars can create an overpressured, expanding bubble, lowering the density in which SNe occur. In addition, expansion of the H ii region will by itself add some momentum to the gas. Neither effect is included in our model, and we would like to understand if this significantly biases our results.

The ionizing luminosity of a cluster of mass M_cluster is Q = 10^49.6M_cluster/(10³ M_⊙) s⁻¹ (Leitherer et al. 1999). These photons will ionize a bubble of gas around the cluster, raising the temperature to 10⁴ K. For density ρ < m_H cm⁻³, this ionized region is not much hotter than the background (which has a temperature only slightly below 10⁴ K, appropriate for warm neutral gas). This means the ionized bubble will not be significantly overpressured compared to the background, so it will not expand significantly. For clusters in backgrounds of density ρ < m_H cm⁻³, we therefore do not expect pre-SN radiative feedback to affect our results.

For higher densities, the equilibrium temperature of the gas is well below 10⁴ K, so the 10⁴ K ionized bubble is significantly overpressured compared to its surroundings. This will allow it to expand, lowering the density in which SNe occur. For uniform density, and neglecting the small range of parameter space where radiation pressure effects will be significant (Krumholz & Matzner 2009), the H ii bubble radius r_II will be governed by the classical Spitzer (1978) solution⁴ (Krumholz 2017, chapter 7),

\begin{equation} r_{\rm II} \approx r_{S,0} \left(\frac{7t}{2\sqrt{3} t_{S,0}}\right)^{4/7}, \end{equation}

(31)

where r_{S, 0} is the Strömgren radius at the start of expansion, given by⁵

\begin{eqnarray} r_{S,0} &=& \left[\frac{3 Q \mu ^2 m_{\rm H}^2}{4(1.1)\pi \alpha _{\rm B} \rho ^2}\right]^{1/3} \nonumber \\ &=& 3.1 Q_{49}^{1/3} n_2^{-2/3} T_{\mathrm{II}, 4}^{0.272 + 0.007 \ln T_{\mathrm{II}, 4}} \mbox{ pc}, \end{eqnarray}

(32)

and where Q₄₉ = Q/10⁴⁹ s⁻¹, μ = 1.33, n₂ = ρ/(μm_H)/100 cm⁻³, T_II is the temperature of the ionized gas, T_{II, 4} = T_II/10⁴ K, |$\alpha _{\rm B} \approx 2.54 \times 10^{-13} T_{\mathrm{II},4}^{-0.8163 - 0.0208 \ln T_{\mathrm{II}, 4}}$| cm³ s⁻¹ is the case B recombination coefficient (Draine 2011), t_{S, 0} = r_{S, 0}/c_II,

\begin{equation} c_{\rm II} = \sqrt{2.2 \frac{k_\mathrm{B} T_{\rm II}}{\mu m_{\rm H}}} = 12 T_{\rm II,4}^{1/2}\mbox{ km s}^{-1}, \end{equation}

(33)

is the ionized gas sound speed. The corresponding density inside the H ii region is

\begin{equation} \rho _{\rm II} = \rho \left(\frac{r}{r_{S,0}}\right)^{-3/2} = \rho \left(\frac{7t}{2\sqrt{3} t_{\mathrm{S},0}}\right)^{-6/7}. \end{equation}

(34)

The mass and velocity of the swept-up shell are

\begin{eqnarray} v_{\rm II} & = & \frac{4}{7} \frac{r_{\rm II}}{t} \end{eqnarray}

(35)

\begin{eqnarray} M_{\rm sh} & = & \frac{4}{3} \pi r_{\rm II}^3 \rho . \end{eqnarray}

(36)

We are interested in the properties of the H ii region at a time of ≈4 Myr, when the first SN occurs; these are

\begin{eqnarray} r_{\rm II} = 27 M_{\rm cluster,3}^{1/7} n_2^{-2/7} \left(\frac{t}{4\mbox{ Myr}}\right)^{4/7} T_{\mathrm{II},4}^{0.402 + 0.003 \ln T_{\mathrm{II}, 4}}\mbox{ pc} \end{eqnarray}

(37)

\begin{eqnarray} \frac{\rho _{\rm II}}{\rho } = 0.077 M_{\rm cluster,3}^{2/7} n_2^{-4/7} \left(\frac{t}{4\mbox{ Myr}}\right)^{-6/7} T_{\mathrm{II},4}^{-0.195 + 0.006\ln T_{\mathrm{II}, 4}} \end{eqnarray}

(38)

\begin{eqnarray} p_{\rm II} & = & 7.7\times 10^5 M_{\rm cluster,3}^{4/7} n_2^{-1/7} \left(\frac{t}{4\mbox{ Myr}}\right)^{-3/7} \nonumber \\ && \times T_{\mathrm{II},4}^{1.41 + 0.02 \ln T_{\mathrm{II}, 4}} \,\mathrm{M}_{\odot }\mbox{ km s}^{-1}, \end{eqnarray}

(39)

where M_{cluster, 3} = M_cluster/10³ M_⊙ and p_II is the momentum of the shell. For the rest of this section, we assume T_{II, 4} = 1.

Based on these results, for M_cluster ≈ 100 M_⊙, the extra momentum in the H ii region (≈2 × 10⁵ M_⊙ km s⁻¹) is ∼50 per cent smaller than that injected by SNe (≈3 × 10⁵ M_⊙ km s⁻¹), and the fractional contribution drops fairly rapidly for more massive clusters thanks to the sublinear scaling of p_II with M_cluster (equation 39). Thus, the extra momentum injected directly by the H ii region is mostly negligible.

The second effect, lowering the density of the medium into which the SNe expand, matters if the region of lowered density encompasses where an SNR would otherwise experience significant cooling: beyond the shell formation radius of the first SN (when the remnant exits the Sedov phase). This radius is (Kim & Ostriker 2015, their equation 8, assuming an energy budget of 10⁵¹ erg per 100 M_⊙ of stars)

\begin{equation} r_{\rm sh} = 6.4 M_{\rm cluster,3}^{0.29} n_2^{-0.42}\mbox{ pc}. \end{equation}

(40)

Given the weak scaling of both r_sh and r_II with M_cluster and n, and the lower coefficient for r_sh, this means that the presence of an H ii region can at least potentially affect the evolution at densities n ≳ 1 cm⁻³.

To quantify the effect of this radiative feedback, we can use the results of Walch & Naab (2015), who ran simulations with and without pre-SN ionization for N_SNe = 1 and ρ ≈ 100 m_H cm⁻³. They found that ionization led to a pre-SN bubble of density ρ ≈ 10 m_H cm⁻³, which resulted in a final momentum 50 per cent higher than their simulation without ionization (runs ‘HCI’ and ‘HC’, respectively). This change in momentum can be explained well by rescaling the results from the non-ionized run to the density within the bubble of the ionized run [using the p ∝ ρ^−1/7 scaling of Cioffi et al. (1988)], with the added extra momentum injected by the H ii region directly. Thus, for single SNe, it appears that the relevant density for a momentum feedback model is the ionized bubble density rather than the background density.

Assuming that this conclusion can be extended to the case of multiple SNe, we can quantity the effects of pre-SN ionizing radiation simply by combining the density scaling in equation (38) with our best-fitting density dependences, p ∝ ρ^−0.06 in the few-SN regime and p ∝ ρ^0.14 in the superbubble regime. In the few-SN regime, this implies an increase in the momentum yield per SN by a factor of

\begin{equation} f_{\rm II,few} \approx 1.17 M_{\rm cluster,3}^{-0.017} n_2^{0.034}. \end{equation}

(41)

The corresponding effect in the superbubble regime is a decrease in the momentum yield by a factor of

\begin{equation} f_{\rm II,many} \approx 0.70 M_{\rm cluster,3}^{0.040} n_2^{-0.080}. \end{equation}

(42)

Thus, in general we expect that the effect of a pre-SN H ii region will be to alter the final momentum yield at the tens of per cent level, with the sign of the effect depending on the whether we are in the few-SN or the superbubble regime.

5.3.2 Type Ia SNe

We just discussed pre-SN alterations to our results – what about changes after the core-collapse SNe occur? In particular, could subsequent Type Ia SNe rejuvenate an old superbubble?

To test this, we rerun a subset of our highest and lowest density simulations with Type Ia SNe added after all the core-collapse SNe have occurred. We add 9.75 × 10⁻⁴ Type Ia SNe per M_⊙ of stars, with the exact number sampled from a Poisson distribution; this rate is taken from Kim et al. (2014), rescaled to a Kroupa (2002) IMF. We draw Type Ia SN delay times from a t⁻¹ distribution, beginning at t = 40 Myr and extending to 100 Myr. This is not meant to be a complete accounting of the full effect of Type Ia SNe; most Type Ia SNe occur after much longer delays, and within 40–100 Myr the delay time distribution is poorly constrained (Maoz, Mannucci & Brandt 2012); this is simply to test the effects of SNe which might occur while the bubble still exists.

We show results for the low- and high-density runs in Figs 18 and 19. For relatively short delay Type Ia SNe, we find that they are consistent with our model if N_SNe is increased accordingly (N_SNe = N_{core-collapse} + N_TypeIa). For long-delay SNe, we caution against using our model; it is not guaranteed that both the progenitor will remain within the cluster and that the bubble will survive much longer than 100 Myr.

Figure 18.

Comparison of simulations with just core-collapse SNe (marked by a blue +), and simulations with core-collapse and Type Ia SNe (marked by a red ×) for a set of simulations with Z = Z_⊙ and ρ = 1.33 × 10⁻³m_H cm⁻³.

Open in new tab Download slide

Same as Fig. 18, except now for ρ = 1.33 × 102mH cm−3 clusters.

Figure 19.

Same as Fig. 18, except now for ρ = 1.33 × 10²m_H cm⁻³ clusters.

Open in new tab Download slide

5.3.3 Self-gravity

Previous studies of SNRs and superbubbles have differed in regard to whether they include or exclude gravitational forces. For example, Martizzi et al. (2015) and Kim & Ostriker (2015) do not include gravity, while Thornton et al. (1998) and Walch & Naab (2015) do include self-gravity. We chose not to include any gravitational forces in our main simulations, and now we estimate what effect that has on our results. In this section, we focus on the effects of the self-gravity of the simulated gas, rather than external gravitational forces.

First, it is useful to understand why we did not include gravity in our simulations. This is partly a philosophical choice. The momentum budget we are attempting to compute is frequently used as an input to models of feedback-regulated star formation or wind generation. In such models, the feedback is compared to the force of gravity either analytically or numerically. For this purpose, the momentum budget in which we are interested is that before the effects of gravity are applied; to include gravity would be in effect to double-count it, by inserting it once into the subgrid feedback model and then a second time into the overall model. However, there is also a practical reason that we omit gravity. Strictly speaking, the self-gravitational potential of an infinite, uniform medium is undefined. We could assume that an external potential dominates, but only if it were spherically symmetric about the cluster, which rules out many potentials of interest, like a galactic potential. Additionally, including a gravitational force, especially self-gravity, would cause our uniform background to collapse. This was not a problem for Thornton et al. (1998), who simulated much less than a free-fall time, but caused Walch & Naab (2015) to limit their simulations to 1 Myr in duration. Unfortunately, limiting the duration of our simulations was not an option; we needed to simulate SNe over a period of at least 30 Myr, much longer than a free-fall time for most of our densities. Rather than artificially require a pressure gradient which ensures equilibrium, we chose to exclude gravitational forces.

Still, gravity exists in real systems, so we should try to understand its effect on our work. First, we will use analytic, simplified arguments to predict the effects of self-gravity on our simulations. We then use a simplified prescription to include self-gravity directly in our numeric simulations. By comparing those results, we can begin to understand the effects of self-gravity and the limitations of our analytic model.

For an arbitrarily thin shell dominated by mass swept up from a constant density background, the force of self-gravity is

\begin{equation} F_\mathrm{grav} = \frac{G M_\mathrm{shell}^2}{2 R_\mathrm{shell}^2} = \frac{8 \pi ^2}{9} G \rho _0^2 R_\mathrm{shell}^4. \end{equation}

(43)

When this inward force becomes greater than the force exerted by the pressure of the hot bubble, the momentum will stop increasing, ending the simulation (unless the bubble would have already mixed into the ISM and the simulation is already completed). This competing force from the hot bubble gas is

\begin{equation} F_\mathrm{gas} \approx 4 \pi R_\mathrm{shell}^2 P_\mathrm{bubble} = 3 (\gamma -1) \frac{E_\mathrm{R,int}}{R_\mathrm{shell}}. \end{equation}

(44)

These two forces become equal at a radius

\begin{eqnarray} R_\mathrm{max} = \left( \frac{27}{8 \pi ^2} (\gamma -1) \frac{ E_\mathrm{R,int}}{G \rho _0^2} \right)^{1/5} \end{eqnarray}

(45)

\begin{eqnarray} \qquad\approx 1200 \left(\frac{N_\mathrm{SNe}}{1000} \right)^{1/5} \left(\frac{\rho _0}{1.33 \times m_\mathrm{H} \,{\rm cm}^{-3}} \right)^{-2/5} {\rm pc,} \end{eqnarray}

(46)

where R_max is the maximum radius the bubble could reach, having assumed that all the injected energy is retained as internal energy (E_{R, int} = N_SNe × 10⁵¹ erg). This provides a simple way to include the effects of gravity via post-processing: using the R_max determined by equation (45), we can truncate a simulation at that radius using the data in Table 2. If the bubble mixes into the ISM at a radius smaller than R_max, then gravity is assumed to have no effect.

We can also rerun a subset of simulations including gravity explicitly. To do this, we calculate the force of self-gravity and the force due to a central point source of mass M_cluster, and use these forces to compute the appropriate source terms for momentum and energy. These source terms are then only applied to the shocked gas (cells with r < R_shock, where R_shock is the radius of the overdensity furthest from the centre). By only applying gravity to shocked gas, we are able to approximate our analytic approach (which only considers gravity of the shell), and avoid the problem of our background collapsing. This effectively assumes that the background is kept in equilibrium by a corresponding thermal or dynamic pressure gradient. We could have explicitly included this pressure gradient in our simulations, but chose not to, so that our with-gravity simulations would better correspond to our without-gravity simulations. A more sophisticated treatment of gravity using higher dimensional simulations would be worthwhile.

In Fig. 20, we compare the results for simulations with and without gravity, and the results of our post-processing model which truncates the without-gravity simulations. We see that including gravity in simulations can decrease the final momentum, and that this decrease can be significant (compared to the uncertainties in our without-gravity model), and as large as a factor of 2.3. We also see that the post-processing truncation model typically overpredicts the final momentum (underpredicting the effects of gravity), relative to the simulations which incorporate gravity directly. This is expected; our post-processing model assumes that all the SN energy remains as thermal energy, but in our simulations some SN energy is converted into kinetic energy, some into potential energy, and most is radiated away. This causes simulations to stop at smaller radii than predicted, leading to lower final momenta than the post-processing model predicts.

Figure 20.

Comparison of simulations with no gravity (marked by a blue +) and simulations with gravity applied to the remnant (marked by a red ×) for all of our simulations with Z = Z_⊙ and ρ = 1.33 × 10⁻³m_H cm⁻³. We also include a simple post-processing prediction for the effect of self-gravity, which truncates our without-gravity simulations if and when they reach R_max (equation 45).

Open in new tab Download slide

We find that self-gravity can significantly change the final momentum, especially for clusters of many SNe, but that momentum is nevertheless still enhanced by a factor of about 4 compared to the single SN case.

5.3.4 Galactic environment

In Section 5.3.3, we investigated the effects of self-gravity, but in some cases a galactic gravitational potential might play a larger role in shaping the late-time dynamics of large bubbles. In this section, we will estimate the effects of a galactic gravitational potential, as well as the effects of rotational shear and disc breakout. For each of these effects, we will use a post-processing correction similar to that used in Section 5.3.3: calculate a radius or time where our assumptions break down, and then use the data in Table 2 to truncate the bubble evolution at that radius or time.

First, we consider the gravitational force perpendicular to a galactic disc, produced by the galactic gravitational potential. The acceleration due to this force can be written as

\begin{equation} g = \frac{z}{r_\mathrm{g}} \frac{v_\mathrm{c}^2}{r_\mathrm{g}}, \end{equation}

(47)

where v_c is the circular velocity at a galactocentric radius r_g, and z is the distance from the mid-plane. For simplicity, we will assume that all of the mass of the shell is in a plane at height z = R_shock. This results in a force

\begin{equation} F_\mathrm{grav} = M_\mathrm{shell} \frac{R_\mathrm{shell}}{r_\mathrm{g}} \frac{v_\mathrm{c}^2}{r_\mathrm{g}}. \end{equation}

(48)

For a Milky Way-like galaxy with v_c = 200 km s⁻¹ and a cluster at r_g = 3 kpc, this gravitational force is equal to the force from the bubble pressure when the shock is at a radius

\begin{eqnarray} R_\mathrm{max} = \left( \frac{9 (\gamma -1)}{4 \pi } \frac{1}{\rho _0} \frac{v_\mathrm{c}^2}{r_\mathrm{g}^2} E_\mathrm{R, int} \right)^{1/5} \end{eqnarray}

(49)

\begin{eqnarray} &&\,\,\,\quad\approx 700 \left( \frac{N_\mathrm{SNe}}{1000} \right)^{1/5} \left( \frac{\rho _0}{1.33 \times m_\mathrm{H}\, {\rm cm}^{-3}} \right)^{1/5} \nonumber \\ && \quad \quad \times \left( \frac{v_\mathrm{c}}{200\, {\rm km\, s}^{-1}} \right)^{-2/5} \left( \frac{r_\mathrm{g}}{3\, {\rm kpc}} \right)^{2/5} {\rm pc}. \end{eqnarray}

(50)

As with our self-gravity model, we create a model which truncates our simulations if and when they reach this radius. In Fig. 21, we compare the results of this galactic gravity model to a subset of our simulations. Similar to the self-gravity results shown in Fig. 20, galactic gravity has the largest effect for the largest clusters. We also find that this galactic gravity model is able to reduce the momentum by a greater factor than our self-gravity model; the self-gravity model never reduced the momentum by more than a factor of 1.4, whereas the galactic gravity model can reduce the momentum by up to a factor of 3.3. As with self-gravity, however, we caution that for many applications, the momentum budget of interest will be that excluding the effects of disc gravity, since if disc gravity is included in the overall model, it should not be double-counted by also including it in the subgrid feedback model.

Figure 21.

Comparison of simulations with no galactic gravity (marked by a blue +) and our simple post-processing model of galactic gravity in a disc for all of our simulations with Z = Z_⊙ and ρ = 1.33 × 10⁻³m_H cm⁻³.

Open in new tab Download slide

Our disc gravity model is incomplete; among other things, it ignores the density and pressure gradients which result from this gravitational potential. As a superbubble expands vertically, it will find less dense and lower pressure gas; as it expands horizontally, it will not experience such gradients in the background gas. This will break the spherical symmetry of the bubble expansion, and can lead to disc breakout where the bubble expansion becomes predominantly vertical. When studying this phenomenon with hydrodynamic simulations, Mac Low, McCray & Norman (1989) found that spherical symmetry is often broken and shell mixing instabilities are significant by the time the bubble has expanded 3–4 scaleheights in the vertical direction. Applying this to a Milky Way-like galaxy, using a thin-disc scaleheight of 100 pc, we can create a model which cuts off the momentum growth when a bubble reaches R_max = 400 pc. The results of this model can be seen in Fig. 22. As with our previous models, the effect of this model is stronger for larger clusters, but this model predicts effects which increase much more rapidly as cluster size increases. While there is no significant effect on our 11-SN simulation, it has the largest effect for the 1008-SN simulation for all of the models we have tested, decreasing the momentum by a factor of about 10. It is important to understand that this behaviour is not a result of the total momentum decreasing for large bubbles. It is simply that, for the largest clusters, the bubble expands to the breakout radius before a significant fraction of the SNe occur, and, by assumption, these additional SNe then contribute no additional momentum. This causes the average momentum per SN to fall, because the extra SNe are counted in the denominator but not the numerator of our average.

Figure 22.

Comparison of simulations with no disc breakout (marked by a blue +) and our simple post-processing model of disc breakout for all of our simulations with Z = Z_⊙ and ρ = 1.33 × 10⁻³m_H cm⁻³.

Open in new tab Download slide

In addition to considering distortions caused by expansion perpendicular to the disc plane, we can also consider distortions due to shear within the disc plane. Adopting a shear time-scale

\begin{equation} t_\mathrm{shear} = \Omega _\mathrm{orb}^{-1} \frac{r_\mathrm{g}}{R_\mathrm{shell}} \end{equation}

(51)

for a disc orbital frequency Ω_orb = v_c/r_g, we can truncate our bubble evolution when the bubble age (t − t_firstSN) exceeds the shear time. Results can be seen in Fig. 23. Once again, only the largest clusters are significantly affected, with the largest effect being a decrease of a factor of 1.75 for the 1008-SN cluster.

Figure 23.

Comparison of simulations with no rotational shear (marked by a blue +) and our simple post-processing model of rotational shear in a disc for all of our simulations with Z = Z_⊙ and ρ = 1.33 × 10⁻³m_H cm⁻³.

Open in new tab Download slide

Comparing these post-processing models, we find that all of the galactic models (gravity, breakout, and shear) are predicted to have larger effects than the self-gravity model, but these galactic models introduce a number of free parameters which we have not explored (disc circular velocity, galactocentric radius, and disc scaleheight). It is important to remember that these models are only very rough estimates. When investigating self-gravity, we found that directly including self-gravity in simulations led to a much stronger effect than predicted by our simple model. The same is likely true for these galactic effects, which we cannot directly include in our simulations. Still, even though every model was structured to decrease the final momentum, we always saw that clustering can lead to an increased momentum efficiency, compared to our single SN results. Moreover, all of these corrections are relatively small for the most efficient clusters, those producing ∼10 SNe (M_cluster ∼ 1000 M_⊙).

Overall, this suggests that superbubble models are less cleanly separable from galactic dynamics than single SNR models. Single, isolated SNRs might only expand to 100 pc over 2 Myr, allowing us to largely ignore effects like disc shear and galactic gravity. As we have shown, these effects can play a significant role for models of clustered SNe and superbubbles. Testing these effects self-consistently goes beyond the capability of our 1D simulations, but we have already gained some insight from our simplified post-processing models. By design, all of these models lowered the final momenta of our bubbles (by as much as a factor of 10, for the disc breakout model applied to our largest bubble), but every model still predicted that clustered SNe lead to enhanced momentum feedback (for instance, our disc breakout model still predicts an average momentum per SN four times larger than the single SN value). Both clustering of SNe and effects from the host galaxy seem to play significant roles in the overall SN momentum budget; it would be useful for higher dimensional simulations to explore these two effects simultaneously in the future.

6 CONCLUSIONS

We perform several hundred 1D simulations to study the momentum delivered to the ISM by clustered SN explosions over a wide range of star cluster sizes, gas densities, and metallicities. Our simulations use a realistic IMF paired with realistic stellar lifetimes, and we evolve them at very high numerical resolution until the momentum of the expanding shell reaches a maximum. At the end of our simulations, we find that our clusters typically retain 1–10 per cent of the injected SN energy (similar to isolated SNe), but clustered SNe produce significantly more momentum per SN than isolated SNe, i.e. the momentum injected by a star cluster is a superlinear function of the number of SNe which explode within it. Clustering has the largest impact for 10–100 SNe, leading to an order-of-magnitude increase in the momentum per SN. When integrating over the observed cluster mass function, our findings suggest that the mean SN momentum budget per mass of stars formed is p/m_* ≈ 1–2 × 10⁴ km s⁻¹, which is ∼0.5–1 dex higher than the value of ≈3000 km s⁻¹ most commonly adopted in the literature.

The increased momentum budget for SNe will have strong implications for any simulation or analytic model which relies on SN momentum injection. In galaxies where the overall star formation rate is high enough for clustering of SNe to be significant for quenching star formation and producing galactic winds (i.e. perhaps not in dwarfs, but almost certainly in more massive galaxies), using our updated SN momentum budget may cause these models’ predicted star formation rates to decrease by the same factor of ∼0.5–1 dex by which the SN momentum budget increases. This may render the models inconsistent with the observed relationship between gas and star formation rate, which will require changes in other free parameters to bring the models back into agreement. The implications for galactic wind launching are less clear, since the increased SN momentum budget will be offset by overall lower star formation rates.

To facilitate the implementation of our results in 3D numerical simulations which include explicit SN momentum injection, we provide a fitting formula for the momentum per SN as a function of cluster size, ambient density, and metallicity. This formula is suitable for implementation in galaxy simulations capable of resolving individual star clusters, typically ∼10²–10⁵ M_⊙. We also provide tabulated outputs from our simulations, for those who wish to calibrate subgrid models at a range of size scales. For lower resolution simulations, we recommend the value of p/m_* ≈ 1–2 × 10⁴ km s⁻¹ we obtain by integrating over the cluster mass function. Regardless of resolution, it is clear that the existing body of simulations may need to be revisited, and some of the strong assumptions which previous authors have adopted to make feedback more effective (e.g. assuming efficient trapping of infrared radiation) may be relaxed. Other numerical schemes, such as turning off radiative cooling for an extended period of time (e.g. Stinson et al. 2006) or stochastically injecting thermal energy in order to delay cooling (e.g. Dalla Vecchia & Schaye 2012), may prove to be closer to reality than had previously been assumed.

Properly capturing the effect of clustering in a 1D simulation requires very high numerical resolution to avoid overcooling through numerical mixing. The resolution requirements may be less severe, and the momentum injection rate lower, in higher dimensions where instabilities may produce mixing at large physical scales. However, the size scale and mixing rate due to instabilities likely depend strongly on the properties of the host galaxy, the nature of the background into which the SNR is propagating, and the microphysical details like magnetic fields and conduction near the contact discontinuity. Since we are unable to directly include these effects in our simulations, we create simple models to estimate the strengths of some of these effects. While these models are able to decrease the momentum of the most massive cluster by a factor of 10, lower mass (more common) clusters are less affected. For each of our models, the average momentum per SN remains greater than the fiducial, unclustered value by a factor of 4, assuming a realistic cluster mass distribution. But these are just simple estimates; no present simulation includes all these effects, and thus the correct momentum budget for clustered SNe in multiple dimensions remains uncertain. While the correct momentum budget may not be as high as ∼0.5–1 dex greater than the commonly adopted single SN value (as we find in one dimension), it is still likely greater than the single SN value. There is clearly an urgent need for further study.

Acknowledgments

We thank Justin Brown, John Forbes, Anna Rosen, and Brant Robertson for useful discussions. We also thank the referee, Norm Murray, for a constructive and insightful report. This work was supported by the NSF through grants AST-1405962 (ESG, MRK, and AD), AST-1229745 (PM), and DGE 1339067 (ESG), by the Australian Research Council through grant ARC DP160100695 (MRK) and by NASA through grant NNX12AF87G (PM). AD's work is also supported by the grants ISF 124/12, I-CORE Programme of the PBC/ISF 1829/12, and BSF 2014-273. MRK thanks the Simons Foundation, which contributed to this work through its support for the Simons Symposium ‘Galactic Superwinds: Beyond Phenomenology’. PM thanks the Préfecture of the Ile-de-France Region for the award of a Blaise Pascal International Research Chair, managed by the Fondation de l'Ecole Normale Supérieure.

1

Source code available at: github.com/egentry/clustered_SNe

2

The situation is more complex for models which include regulation of star formation by FUV radiation instead of or in addition to SN feedback (e.g. Krumholz, McKee & Tumlinson 2009; Ostriker, McKee & Leroy 2010; Krumholz 2013). In these models, the effects of enhanced momentum injection will be more modest or absent, depending on the details of the individual model.

3

There exist prescriptions for approximating mixing instabilities in 1D codes (e.g. Duffell 2016), but none were incorporated in our work.

4

This solution assumes t ≫ t_{S, 0}, otherwise expansion is expected to be negligible. For the cluster masses considered, this assumption typically fails for n < 1 cm⁻³, but we already expected minimal expansion in such backgrounds.

5

Note that this and the following expressions assume that He is singly ionized.

REFERENCES

Bryan

G. L.

et al.

2014

,

ApJS

,

211

,

19

Crossref

Search ADS

Castor

J.

,

McCray

R.

,

Weaver

R.

1975

,

ApJ

,

200

,

L107

Crossref

Search ADS

Chevalier

R. A.

1974

,

ApJ

,

188

,

501

Crossref

Search ADS

Cioffi

D. F.

,

McKee

C. F.

,

Bertschinger

E.

1988

,

ApJ

,

334

,

252

Crossref

Search ADS

Cowie

L. L.

,

McKee

C. F.

1977

,

ApJ

,

211

,

135

Crossref

Search ADS

Creasey

P.

,

Theuns

T.

,

Bower

R. G.

2013

,

MNRAS

,

429

,

1922

Crossref

Search ADS

da Silva

R. L.

,

Fumagalli

M.

,

Krumholz

M.

2012

,

ApJ

,

745

,

145

Crossref

Search ADS

da Silva

R. L.

,

Fumagalli

M.

,

Krumholz

M. R.

2014

,

MNRAS

,

444

,

3275

Crossref

Search ADS

Dalla Vecchia

C.

,

Schaye

J.

2012

,

MNRAS

,

426

,

140

Crossref

Search ADS

Dekel

A.

,

Krumholz

M. R.

2013

,

MNRAS

,

432

,

455

Crossref

Search ADS

Dekel

A.

,

Silk

J.

1986

,

ApJ

,

303

,

39

Crossref

Search ADS

Draine

B. T.

2011

,

Physics of the Interstellar and Intergalactic Medium

,

Princeton Univ. Press

,

Princeton, NJ

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Duffell

P. C.

2016

,

ApJ

,

821

,

76

Crossref

Search ADS

Ekström

S.

et al.

2012

,

A&A

,

537

,

A146

Crossref

Search ADS

Ertl

T.

,

Janka

H.-T.

,

Woosley

S. E.

,

Sukhbold

T.

,

Ugliano

M.

2016

,

ApJ

,

818

,

124

Crossref

Search ADS

Faucher-Giguère

C.-A.

,

Quataert

E.

,

Hopkins

P. F.

2013

,

MNRAS

,

433

,

1970

Crossref

Search ADS

Ferland

G. J.

,

Korista

K. T.

,

Verner

D. A.

,

Ferguson

J. W.

,

Kingdon

J. B.

,

Verner

E. M.

1998

,

PASP

,

110

,

761

Crossref

Search ADS

Fierlinger

K. M.

,

Burkert

A.

,

Ntormousi

E.

,

Fierlinger

P.

,

Schartmann

M.

,

Ballone

A.

,

Krause

M. G. H.

,

Diehl

R.

2016

,

MNRAS

,

456

,

710

Crossref

Search ADS

Gerritsen

J. P. E.

1997

,

PhD thesis

,

Groningen University

,

the Netherlands

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Goldbaum

N. J.

,

Krumholz

M. R.

,

Forbes

J. C.

2016

,

ApJ

,

827

,

28

Crossref

Search ADS

Haardt

F.

,

Madau

P.

2012

,

ApJ

,

746

,

125

Crossref

Search ADS

Hayward

C. C.

,

Hopkins

P. F.

2016

,

MNRAS

in press

Hennebelle

P.

,

Iffrig

O.

2014

,

A&A

,

570

,

A81

Crossref

Search ADS

Hopkins

P. F.

,

Quataert

E.

,

Murray

N.

2011

,

MNRAS

,

417

,

950

Crossref

Search ADS

Hopkins

P. F.

,

Quataert

E.

,

Murray

N.

2012

,

MNRAS

,

421

,

3522

Crossref

Search ADS

Hopkins

P. F.

,

Kereš

D.

,

Oñorbe

J.

,

Faucher-Giguère

C.-A.

,

Quataert

E.

,

Murray

N.

,

Bullock

J. S.

2014

,

MNRAS

,

445

,

581

Crossref

Search ADS

Iffrig

O.

,

Hennebelle

P.

2015

,

A&A

,

576

,

A95

Crossref

Search ADS

Jenkins

E. B.

,

Tripp

T. M.

2011

,

ApJ

,

734

,

65

Crossref

Search ADS

Katz

N.

1992

,

ApJ

,

391

,

502

Crossref

Search ADS

Keller

B. W.

,

Wadsley

J.

,

Benincasa

S. M.

,

Couchman

H. M. P.

2014

,

MNRAS

,

442

,

3013

Crossref

Search ADS

Kim

J.-h.

et al.

2014

,

ApJS

,

210

,

14

Crossref

Search ADS

Kim

C.-G.

,

Ostriker

E. C.

2015

,

ApJ

,

802

,

99

Crossref

Search ADS

Kim

C.-G.

,

Kim

W.-T.

,

Ostriker

E. C.

2011

,

ApJ

,

743

,

25

Crossref

Search ADS

Kim

C.-G.

,

Ostriker

E. C.

,

Raileanu

R.

2016

preprint (arXiv:1610.03092)

Kimm

T.

,

Cen

R.

,

Devriendt

J.

,

Dubois

Y.

,

Slyz

A.

2015

,

MNRAS

,

451

,

2900

Crossref

Search ADS

Krause

M.

,

Fierlinger

K.

,

Diehl

R.

,

Burkert

A.

,

Voss

R.

,

Ziegler

U.

2013

,

A&A

,

550

,

A49

Crossref

Search ADS

Kroupa

P.

2002

,

Science

,

295

,

82

Crossref

Search ADS

PubMed

Krumholz

M. R.

2013

,

MNRAS

,

436

,

2747

Crossref

Search ADS

Krumholz

M. R.

2014

,

Phys. Rep.

,

539

,

49

Crossref

Search ADS

Krumholz

M. R.

2017

,

Star Formation

,

World Scientific Press

,

Singapore

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Krumholz

M. R.

,

Matzner

C. D.

2009

,

ApJ

,

703

,

1352

Crossref

Search ADS

Krumholz

M. R.

,

McKee

C. F.

,

Tumlinson

J.

2009

,

ApJ

,

699

,

850

Crossref

Search ADS

Krumholz

M. R.

,

Fumagalli

M.

,

da Silva

R. L.

,

Rendahl

T.

,

Parra

J.

2015

,

MNRAS

,

452

,

1447

Crossref

Search ADS

Leitherer

C.

et al.

1999

,

ApJS

,

123

,

3

Crossref

Search ADS

Mac Low

M.-M.

,

McCray

R.

,

Norman

M. L.

1989

,

ApJ

,

337

,

141

Crossref

Search ADS

Maoz

D.

,

Mannucci

F.

,

Brandt

T. D.

2012

,

MNRAS

,

426

,

3282

Crossref

Search ADS

Markevitch

M.

,

Vikhlinin

A.

2007

,

Phys. Rep.

,

443

,

1

Crossref

Search ADS

Martizzi

D.

,

Faucher-Giguère

C.-A.

,

Quataert

E.

2015

,

MNRAS

,

450

,

504

Crossref

Search ADS

Murray

N.

,

Quataert

E.

,

Thompson

T. A.

2005

,

ApJ

,

618

,

569

Crossref

Search ADS

Noh

W. F.

1987

,

J. Comput. Phys.

,

72

,

78

Crossref

Search ADS

Ostriker

E. C.

,

Shetty

R.

2011

,

ApJ

,

731

,

41

Crossref

Search ADS

Ostriker

E. C.

,

McKee

C. F.

,

Leroy

A. K.

2010

,

ApJ

,

721

,

975

Crossref

Search ADS

Pejcha

O.

,

Thompson

T. A.

2015

,

ApJ

,

801

,

90

Crossref

Search ADS

Sharma

P.

,

Roy

A.

,

Nath

B. B.

,

Shchekinov

Y.

2014

,

MNRAS

,

443

,

3463

Crossref

Search ADS

Shetty

R.

,

Ostriker

E. C.

2012

,

ApJ

,

754

,

2

Crossref

Search ADS

Spitzer

L.

Jr

1978

,

J. R. Astron. Soc. Can.

,

72

,

349

Stinson

G.

,

Seth

A.

,

Katz

N.

,

Wadsley

J.

,

Governato

F.

,

Quinn

T.

2006

,

MNRAS

,

373

,

1074

Crossref

Search ADS

Sukhbold

T.

,

Ertl

T.

,

Woosley

S. E.

,

Brown

J. M.

,

Janka

H.-T.

2016

,

ApJ

,

821

,

38

Crossref

Search ADS

Thompson

T. A.

,

Krumholz

M. R.

2016

,

MNRAS

,

455

,

334

Crossref

Search ADS

Thompson

T. A.

,

Quataert

E.

,

Murray

N.

2005

,

ApJ

,

630

,

167

Crossref

Search ADS

Thornton

K.

,

Gaudlitz

M.

,

Janka

H.-T.

,

Steinmetz

M.

1998

,

ApJ

,

500

,

95

Crossref

Search ADS

Toro

E.

,

Spruce

M.

,

Speares

W.

1994

,

Shock Waves

,

4

,

25

Crossref

Search ADS

Walch

S.

,

Naab

T.

2015

,

MNRAS

,

451

,

2757

Crossref

Search ADS

Woosley

S. E.

,

Heger

A.

2007

,

Phys. Rep.

,

442

,

269

Crossref

Search ADS

Woosley

S. E.

,

Heger

A.

2015

,

ApJ

,

810

,

34

Crossref

Search ADS

Yadav

N.

,

Mukherjee

D.

,

Sharma

P.

,

Nath

B. B.

2016

,

MNRAS

in press

SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

Table 1. Overview of numeric results.

Table 2. Momentum evolution.

Please note: Oxford University Press are not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

APPENDIX A: CODE VERIFICATION

A1 Sedov verification

The mass and energy of each SN are injected into the innermost zone, with all the energy injected as thermal energy. This is not a realistic configuration; at no stage do we expect a uniformly mixed sphere, of the order of 0.1 pc in radius, which is overpressured but not yet expanding. Given these convenient but unphysical initial conditions, we need to verify that our system will evolve into a realistic configuration.

We can look at an early time snapshot of a single SN simulation to verify that the system accurately relaxes into a physical configuration. At early times, cooling losses should still be negligible, so we expect our system to be in the Sedov phase. Fig. A1 shows a snapshot of our numerical results, compared with the analytic Sedov prediction. There is a noticeable overdensity at inner radii, but this is to be expected: the analytic Sedov solution assumes no ejecta mass, whereas our simulation includes ejecta mass. The extra ejecta mass appears as an overdensity at inner radii. Excepting that, our simulation is in good agreement with the Sedov prediction so we consider our injection scheme valid.

Comparison of our numeric results (solid line) against the analytic Sedov solution (dashed line) for a ρ = 1.33 × 10−3mH cm−3, Z = Z⊙, Mcluster = 102 M⊙ (NSNe = 1) cluster, at t = 0.17 Myr.

Figure A1.

Comparison of our numeric results (solid line) against the analytic Sedov solution (dashed line) for a ρ = 1.33 × 10⁻³m_H cm⁻³, Z = Z_⊙, M_cluster = 10² M_⊙ (N_SNe = 1) cluster, at t = 0.17 Myr.

Open in new tab Download slide

A2 Thornton et al. verification

We also verified our code against the results of Thornton et al. (1998), who measured the total energy from single SNe. We ran single SN simulations at the same background conditions as Thornton et al., fixing the SN ejecta mass to be 3 M_⊙ with an ejecta metallicity equal to the background metallicity, and extracting results at the same time as Thornton et al. We compare our simulations to the model provided by Thornton et al. in Figs A2 and A3.

Figure A2.

Verification that our code can reproduce the results of Thornton et al. (1998), for total energy contained within the SNR (E_{R, tot}) at the completion time defined by Thornton et al. (1998).

Open in new tab Download slide

Figure A3.

Same as Fig. A2, but now with total energy as a function of metallicity.

Open in new tab Download slide

We judge that our residuals are comparable to the residuals present in the data of Thornton et al., and we are not surprised that there are discrepancies. We use different initial conditions: our simulation injects all of the SN energy into the innermost zone as thermal energy, while that of Thornton et al. spreads the energy across 150 zones, and adds some of it as kinetic energy. We were able to use different initial conditions because we used different hydrodynamic solvers: Thornton et al. use a finite-difference method which cannot handle the strong shock which occurs by injecting all the energy into one zone, while our finite-volume method is much more robust to these strong shock conditions. Finally, we use a cooling package which differs from the cooling function used by Thornton et al. All these differences lead us to expect the minor discrepancies between our results and the results of Thornton et al.

Download all slides

Month:	Total Views:
November 2016	1
January 2017	15
February 2017	13
March 2017	16
April 2017	11
May 2017	8
June 2017	13
July 2017	3
September 2017	1
October 2017	10
November 2017	20
December 2017	14
January 2018	15
February 2018	5
March 2018	7
April 2018	5
May 2018	3
June 2018	14
July 2018	8
August 2018	13
September 2018	8
October 2018	12
November 2018	18
December 2018	8
January 2019	6
February 2019	9
March 2019	12
April 2019	7
May 2019	23
June 2019	16
July 2019	12
August 2019	14
September 2019	23
October 2019	24
November 2019	9
December 2019	10
January 2020	11
February 2020	14
March 2020	7
April 2020	7
May 2020	9
June 2020	15
July 2020	19
August 2020	6
September 2020	10
October 2020	12
November 2020	14
December 2020	8
January 2021	16
February 2021	10
March 2021	18
April 2021	13
May 2021	13
June 2021	5
July 2021	12
August 2021	11
September 2021	7
October 2021	5
November 2021	18
January 2022	6
February 2022	4
March 2022	13
April 2022	7
May 2022	22
June 2022	5
July 2022	13
August 2022	12
September 2022	14
October 2022	9
November 2022	8
December 2022	11
January 2023	2
February 2023	3
March 2023	8
April 2023	4
May 2023	5
June 2023	2
July 2023	10
August 2023	8
September 2023	9
October 2023	9
November 2023	9
December 2023	20
January 2024	12
February 2024	20
March 2024	16
April 2024	14
May 2024	15
June 2024	17
July 2024	24
August 2024	12
September 2024	17
October 2024	14
November 2024	16
December 2024	18
January 2025	18
February 2025	15
March 2025	18
April 2025	6
May 2025	6

Article Contents

Enhanced momentum feedback from clustered supernovae Free

Abstract

1 INTRODUCTION

2 NUMERICAL METHODS

2.1 Initial conditions

2.2 Hydrodynamics

2.3 Cooling

2.4 Cluster model

2.5 SN injection model

2.5.1 Winds

3 NUMERICAL RESULTS

4 THE MOMENTUM BUDGET OF CLUSTERED SNe

4.1 Qualitative analysis

4.1.1 The small-N regime

4.1.2 The superbubble regime

4.2 Quantitative model

5 DISCUSSION

5.1 Implications of high momentum efficiency of clustered SNe

5.2 Comparison to previous work and convergence study

5.3 The effects of additional physics

5.3.1 Pre-SN radiative feedback

5.3.2 Type Ia SNe

5.3.3 Self-gravity

5.3.4 Galactic environment

6 CONCLUSIONS

Acknowledgments

REFERENCES

SUPPORTING INFORMATION

APPENDIX A: CODE VERIFICATION

A1 Sedov verification

A2 Thornton et al. verification

Supplementary data

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only

Enhanced momentum feedback from clustered supernovae