Bursting Bubbles: Clustered Supernova Feedback in Local and High-redshift Galaxies

In O21, we considered a simplified slab geometry for the ISM with a mean mass density of ${\bar{\rho }}_{g}={{\rm{\Sigma }}}_{g}/2H$, where a star cluster forms at the galactic midplane with a mass ${M}_{\mathrm{cl}}=\pi {H}^{2}{{\rm{\Sigma }}}_{g}^{2}/{{\rm{\Sigma }}}_{\mathrm{crit}}$ (according to Grudić et al. 2018, Σ_crit = 2800 M_⊙ pc⁻²). To model the evolution of the superbubble, we take the bubble to be in a momentum-conserving phase, with the momentum of the shock front at a radius R_b , having swept up the mass of gas within that radius, to be ${P}_{b}=\tfrac{4}{3}\pi {R}_{b}^{3}{\bar{\rho }}_{g}{\rm{d}}{R}_{b}/{\rm{d}}t$ (Equation (3), O21, also Fielding et al. 2018; El-Badry et al. 2019), where dR_b /dt ≡ v_b is the expansion velocity of the superbubble. At all points in time, we balance this momentum with the cumulative momentum injected up until that time by SNe from the central star cluster (i. e., P_b = P_SNe). To model the rate and nominal momentum injection of SNe, we invoke a power-law delay time distribution, with ${\rm{d}}{N}_{\mathrm{SN}}/$dt ∝ t^−α (see Appendix A of Orr et al. 2019 and O21 for a more detailed discussion; in this Letter, we take α = 0.46) and assume that a fiducial momentum of ${\left(P/{m}_{\star }\right)}_{0}$ is injected by each SN (≈3000 km s⁻¹, Martizzi et al. 2015, normalized as 1 SN per 100 M_⊙). Notionally, all SNe occur over a time period 0 < t < t_SNe corresponding to the time from first SNe to occur in the star cluster (lifetime of the most massive star formed) until the time of the last SNe to occur (lifetime of the least massive star to undergo a core-collapse SN, ≈ 40 Myr), which divides up our parameter space into cases where the bubble evolves with and without additional SN momentum being injected.

Balancing the bubble and feedback momenta, we found relations for the bubble radius and shock front velocity in time (Equations (7) and (8), O21):

$$\begin{eqnarray*}&&\displaystyle \frac{{R}_{b}}{H}={\left[6\displaystyle \frac{{{\rm{\Sigma }}}_{g}}{{{\rm{\Sigma }}}_{\mathrm{crit}}}\displaystyle \frac{{\left(P/{m}_{\star }\right)}_{0}}{H/{t}_{\mathrm{SNe}}}\right]}^{\tfrac{1}{4}}\left\{\begin{array}{l}\tfrac{1}{{\left(2-\alpha \right)}^{1/4}}{\left(\tfrac{t}{{t}_{\mathrm{SNe}}}\right)}^{\tfrac{2-\alpha }{4}}\\ \qquad \qquad 0\lt t\lt {t}_{\mathrm{SNe}}\\ {\left[\tfrac{t}{{t}_{\mathrm{SNe}}}-\left(\tfrac{1-\alpha }{2-\alpha }\right)\right]}^{\tfrac{1}{4}}\\ \qquad \qquad t\gt {t}_{\mathrm{SNe}}\end{array}\right.\end{eqnarray*}$$

$$\begin{eqnarray*}&&\begin{array}{l}\displaystyle \frac{{v}_{b}}{\sigma }=\displaystyle \frac{1}{4{\rm{\Omega }}{t}_{\mathrm{SNe}}}{\left[6\displaystyle \frac{{{\rm{\Sigma }}}_{g}}{{{\rm{\Sigma }}}_{\mathrm{crit}}}\displaystyle \frac{{\left(P/{m}_{\star }\right)}_{0}}{H/{t}_{\mathrm{SNe}}}\right]}^{\tfrac{1}{4}}\\ \,\ \ \times \,\left\{\begin{array}{l}{\left(2-\alpha \right)}^{\tfrac{3}{4}}{\left(\tfrac{t}{{t}_{\mathrm{SNe}}}\right)}^{-\tfrac{(2+\alpha )}{4}}\\ \qquad \qquad 0\lt t\lt {t}_{\mathrm{SNe}}\\ {\left[\tfrac{t}{{t}_{\mathrm{SNe}}}-\left(\tfrac{1-\alpha }{2-\alpha }\right)\right]}^{-\tfrac{3}{4}}\\ \qquad \qquad t\gt {t}_{\mathrm{SNe}}\end{array}\right.\end{array}\end{eqnarray*}$$

The outcomes of the evolution of these superbubbles are broken down into four cases, depending on whether or not the bubbles break out of the disk (at a time t_BO), and if the central star cluster is still producing SNe. If the superbubble comes into pressure equilibrium with the surrounding ISM (for a turbulent ISM, ${P}_{\mathrm{ISM}}\sim {\bar{\rho }}_{g}{\sigma }^{2}$), it will not maintain coherence in its expansion to reach the gas-disk scale height, and instead the shock front will fragment/stall (as shown in the simulations by Fielding et al. 2018). Thus, we take the superbubbles to stall and fragment when v_b ≈ σ, demarcating the difference between cases where the remnant successfully and unsuccessfully reaches the disk scale height.

We considered the following cases in O21,

• PBO Case: "Powered Breakout," SNe remnant superbubble reaches the gas-disk scale height, R_b = H, before the central star cluster ceases producing SNe, t_BO < t_SNe.
• CBO Case: "Coasting (unpowered) Breakout," remnant reaches the gas disk scale height, R_b = H, after the central star cluster ceases producing SNe, t_BO > t_SNe.
• CF Case: "Coasting (unpowered) Fragmentation," the remnant fragments in the turbulent ISM (i.e., the velocity of the shock front falls below the turbulent velocity of the ISM), v_b < σ, before reaching the gas disk scale height, R_b (v_b = σ) < H, after the central star cluster ceases producing SNe, t(v_b = σ) > t_SNe.
• PS Case: "Powered Stall," bubble expansion stalls in the turbulent ISM, v_b < σ, before reaching the gas-disk scale height, R_b (v_b = σ) < H, before the central star cluster ceases producing SNe, t(v_b = σ) < t_SNe.

Following our assumptions regarding star cluster formation efficiency, and that locally, gas in galaxies finds itself marginally stable against gravitational fragmentation and collapse with Toomre-${\tilde{Q}}_{\mathrm{gas}}\approx 1$, we found in O21 that the boundaries between these four cases could be entirely expressed in terms of local gas fraction ${\tilde{f}}_{g}\equiv {{\rm{\Sigma }}}_{g}/({{\rm{\Sigma }}}_{g}+{{\rm{\Sigma }}}_{\star })$ and inverse dynamical time Ω ≡ v_c /R. We refer to Table 1 for the boundaries between cases, in ${\tilde{f}}_{g}$–Ω space, along with the constraints in Ω derived in O21 (and referencing Equation numbers therein).

Table 1. Superbubble Outcome Boundaries

Boundary Cases (Description)	Boundary Equation	Parametric Constraint	O21 Equation No.
PBO/PS Case ("Powered Breakout/Stall")	${\tilde{f}}_{g}=\tfrac{\sqrt{2}\pi G}{3}\tfrac{{{\rm{\Sigma }}}_{\mathrm{crit}}}{{\left(P/{m}_{\star }\right)}_{0}}{\left(\tfrac{4{\rm{\Omega }}{t}_{\mathrm{SNe}}}{2-\alpha }\right)}^{1-\alpha }\tfrac{1}{{\rm{\Omega }}}$	Ω ≥ (2 − α)/4tSNe	Equation (9)
PBO/CBO Case ("Powered/Coasting Breakout")	${\tilde{f}}_{g}=\tfrac{\sqrt{2}\pi G}{3}\tfrac{{{\rm{\Sigma }}}_{\mathrm{crit}}}{{\left(P/{m}_{\star }\right)}_{0}}\left(\tfrac{2-\alpha }{4{\rm{\Omega }}{t}_{\mathrm{SNe}}}\right)\tfrac{1}{{\rm{\Omega }}}$		Equation (12)
CBO/CF Case("Coasting Breakout/Fragmentation")	${\tilde{f}}_{g}=\tfrac{\sqrt{2}\pi G}{3}\tfrac{{{\rm{\Sigma }}}_{\mathrm{crit}}}{{\left(P/{m}_{\star }\right)}_{0}}\tfrac{1}{{\rm{\Omega }}}$	Ω ≤ (2 − α)/4tSNe	Equation (15)
CF/PS Case("Coasting Fragmentation/Powered Stall")	${\tilde{f}}_{g}=\tfrac{\sqrt{2}\pi G}{3}\tfrac{{{\rm{\Sigma }}}_{\mathrm{crit}}}{{\left(P/{m}_{\star }\right)}_{0}}{\left(\tfrac{4{\rm{\Omega }}{t}_{\mathrm{SNe}}}{2-\alpha }\right)}^{3}\tfrac{1}{{\rm{\Omega }}}$		Equation (13)

Download table as:  ASCIITypeset image

Here we predict the local turbulence driving scale from the fragmentation of superbubbles within the ISM and effective strength of feedback in regions that host bubble breakout for observed galaxies. In O21, we calculated the fragmentation/stall radius for bubbles in the CF/PS case and found (Equations (17) and (18), O21),

$$\begin{eqnarray*}\displaystyle \frac{{R}_{b}({t}_{\mathrm{frag}})}{H}=\left\{\begin{array}{ll}{\left[\tfrac{{\left(P/{m}_{\star }\right)}_{0}}{{{\rm{\Sigma }}}_{\mathrm{crit}}}\tfrac{3{\rm{\Omega }}{\tilde{f}}_{g}}{\sqrt{2}\pi G}\right]}^{\tfrac{1}{3}} & {\bf{CF\; Case}}\\ {\left[\tfrac{{\left(P/{m}_{\star }\right)}_{0}}{{{\rm{\Sigma }}}_{\mathrm{crit}}}\tfrac{3{\rm{\Omega }}{\tilde{f}}_{g}}{\sqrt{2}\pi G}\right]}^{\tfrac{1}{2+\alpha }} & \\ \times {\left(\tfrac{2-\alpha }{4{\rm{\Omega }}{t}_{\mathrm{SNe}}}\right)}^{\tfrac{1-\alpha }{2+\alpha }} & {\bf{PS\; Case}}\end{array}\right.\end{eqnarray*}$$

Similarly, we predicted the effective strength of feedback, ${(P/{M}_{\star })}_{\mathrm{eff}}$, i.e., the fraction of momentum deposited into the ISM versus lost to outflows in the event that the superbubble were to break out of the ISM (PBO case), and found (Equation (20), O21),

$$\begin{eqnarray*}&&\displaystyle \frac{{\left(P/{m}_{\star }\right)}_{\mathrm{eff}}}{{\left(P/{m}_{\star }\right)}_{0}}={\left[\displaystyle \frac{\sqrt{2}\pi G}{3}\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{crit}}}{{\left(P/{m}_{\star }\right)}_{0}}\displaystyle \frac{2-\alpha }{4{\rm{\Omega }}{t}_{\mathrm{SNe}}}\displaystyle \frac{1}{{\tilde{f}}_{g}{\rm{\Omega }}}\right]}^{\left(\tfrac{1-\alpha }{2-\alpha }\right)}.\end{eqnarray*}$$

This necessarily would affect the slope of the Kennicutt–Schmidt (KS) relation Kennicutt & Evans (2012), in a feedback-regulated framework (e.g., Faucher-Giguere et al. 2013), for ISM patches in PBO parameter space.

Figure 2 shows the ${\tilde{f}}_{g}$–Ω parameter space of superbubble outcomes colored by the predicted fragmentation radius and effective strength of feedback, for their appropriate cases, with the same observations as Figure 1. We predict that for local star-forming galaxies (and the conditions of the solar circle), most superbubbles which do not break out of the ISM nevertheless still grow to an appreciable fraction of the gas scale height (R_b (t_frag)/H ≳ 0.7). Indeed this might be an expected attractor state, as if the gas scale height is to be set by turbulence, and the vertical turbulent field is to be driven by SNe, then we ought to expect that SNe have a turbulence driving scale of roughly the scale height.

Zoom In Zoom Out Reset image size

We see that interestingly, the effective strength of feedback is perhaps not dramatically reduced in z ∼ 2 galaxies, but for local supernova-driven outflows, we predict that, e.g., NGC 253 might have only ∼10% of the feedback momentum from the central starburst deposited into its ISM. This suggests that although outflows might be ubiquitous at cosmic noon, their effects are significantly different in regards to the ability to locally regulate the ISM.

Comparing to the simulations of Fielding et al. (2018) that fall in the PS and PBO cases, respectively, we also find satisfactory agreement with our predictions. Their simulated superbubble that failed to break out grew to ∼0.9H, before fragmenting and churning with a relative size of ∼0.7–0.8H, and the simulation that successfully broke out did so after approximately ∼2–3 Myr. This fragmentation scale was slightly larger than we predict here, but their numerical setup slightly differed from our model assumptions, having a vertically stratified inhomogeneous ISM and flat (α = 0) SNe time distribution, which may account for the difference. In the case of the successful breakout simulation, their flat SNe time distribution and t_SNe ≈ 30 Myr may account for the difference between our predicted ${\left(P/{m}_{\star }\right)}_{\mathrm{eff}}$ and the simulation (see Equation (19), O21).

In Figure 3, we plot observed H i bubble radii and expansion velocities from local star-forming galaxies in the things survey by Bagetakos et al. (2011, we used things rotation curves from de Blok et al. 2008 to estimate local H) to interpret the likely outcome of these bubbles. Bagetakos et al. (2011) identified gaps and voids in spatially and velocity-resolved data of the things H i disks, fitting ellipsoids to find H i bubble sizes and expansion velocities. Their sample was divided up into three types of H i bubble; here we consider only their type 3, where the bubble is still intact with both near and far edges detected, comparable with still-evolving superbubbles in our model.

Zoom In Zoom Out Reset image size

Predicting whether or not we expect an observed superbubble to fragment or breakout is possible when considering the t < t_SNe cases of Equations (7) and (8) of O21 (see Section 2.1). Taking the ratio v_b /R_b , we find that this evolves as ${v}_{b}\propto {R}_{b}^{-(2+\alpha )/(2-\alpha )}$. And so, bubbles observed to be below a line of this constant slope divide the v_b –R_b space into bubbles that we expect to fragment and those that we expect to break out.

The observed bubbles whose radii appear to already be greater than H all come from short dynamical time regions, where 1/Ω = t_dyn < 20 Myr = t_SNe/2. This suggests that the H i bubbles are remaining fairly coherent after they have already broken out of the galactic nuclei, or that we are systematically underestimating H in these regions. A number of these bubbles are in the inner ring of NGC 4736, which appears to be a dynamically induced starburst (Munoz-Tunon et al. 2004), for which our assumption of ${\tilde{Q}}_{\mathrm{gas}}\approx 1$ may not hold. As well, identifying intact H i bubbles here may be problematic given the predominantly molecular nature of the central regions of local L_⋆ star-forming galaxies (Jiménez-Donaire et al. 2019).

Considering the observed intact bubbles from regions with Ω < 50 Gyr⁻¹ (for v_c ≈ 200 km s⁻¹, this is R > v_c /Ω ≈ 4 kpc), the ensemble of radii and expansion velocities suggests that all of these H i bubbles will eventually fragment in the ISM rather than drive significant outflows/fountains. This is consistent with the ${\tilde{f}}_{\mathrm{gas}}$-Ω profiles of local star-forming galaxies (see Figure 1), where the only regions that host superbubble breakout are the central starbursts (where various gas dynamics have fed the formation of central super star clusters).

As discussed in Section 3.4 of O21, the primary difficulty in realizing coasting (CBO and CF cases) outcomes appears to lie in the fact that the star-forming extent of the vast majority of galaxies does not reach so far out (radially) to have dynamical times exceeding a 100 Myr (i.e., 1/t_dyn = Ω = 10 Gyr⁻¹). For the most part, rotation curves in galaxies are sufficiently high, and their star-forming edges sufficiently close, that dynamical times remain shorter than 4t_SNe/(2 − α), and only powered outcomes are seen. The only exceptions may be in ultra-diffuse dwarfs (UDGs), having low v_c and large extents (Beasley et al. 2016). Even then, this model is only relevant for those that still maintain some star-forming gas (improbable for UDGs).

Alternative models for star cluster formation efficiency have been proposed, arguing that star formation proceeds at a constant efficiency of roughly 1% per freefall time (Krumholz & Tan 2007; see Krumholz et al. 2019 for a review of (low) star formation efficiency in clusters). If we were to adopt a constant efficiency per Toomre-mass model, as opposed to the Grudić et al. (2018) model, where M_cl = _cl M_g ≈ _cl π H²Σ_g , where _cl = 0.01, and holding the rest of the model fixed, then throughout the text the only difference required would be to replace Σ_crit → Σ_g /_cl. Whereupon, we would plot all our figures in ${\tilde{f}}_{g}/{{\rm{\Sigma }}}_{g}$–Ω space, rather than ${\tilde{f}}_{g}$–Ω space. The rest of the results and analysis would remain unchanged. Figure 4 shows the case boundaries in the somewhat unusual ${\tilde{f}}_{g}/{{\rm{\Sigma }}}_{g}$–Ω space.

Zoom In Zoom Out Reset image size

The main difference of this alternative model is that breakout is now dependent on the local disk surface density, as ${{\rm{\Sigma }}}_{\mathrm{disk}}\approx {{\rm{\Sigma }}}_{\star }+{{\rm{\Sigma }}}_{g}={{\rm{\Sigma }}}_{g}/{\tilde{f}}_{g}$. In fact, such a model would imply that breakout only occurs at lower (relatively speaking) disk surface densities for a given local dynamical time: above a local disk surface density, the resulting superbubble is effectively smothered by the disk. At the critical dynamical time (where all four case boundaries intersect) this local disk surface density is ≈45 M_⊙ pc⁻², and it grows to ≈130 M_⊙ pc⁻² at Ω = 10² Gyr⁻¹. Given that regions with shorter dynamical times are generally coreward in disk galaxies, and that generally the disk (stellar + gaseous) surface densities of the central regions in Milky Way mass galaxies greatly exceed 100 M_⊙ pc⁻² (see the PHANGS-ALMA sample of Sun et al. 2020), this model suggests that superbubble breakout, and thus outflows, would not occur in the inner disk regions of galaxies, only in their outskirts. This is contrary to many observations of galactic winds and outflows, specifically those of galaxies without AGN, which nonetheless tend to report outflows and fountains originating from the inner regions of galaxies (e.g., Bolatto et al. 2013). Consequently, these data disfavor a constant (per Toomre-mass) star cluster formation efficiency within this superbubble feedback model.

In this Letter, we compared a model (derived in O21) of clustered SNe feedback in disk environments with observations of local and high-redshift star-forming galaxies. Of specific interest, we tested our predictions from O21 of the ability of supernova-driven superbubbles to break out of the gas disk of a galaxy with known hosts of superwinds, and galaxies thought to lack them. As well, we examined observed H i bubbles/holes in the context of our predicted scalings for bubble radii and velocities.

Key takeaways from comparing this model to observations include:

• 1.  
  Spatially resolved observations of z ≈ 0 star-forming galaxies suggest that most star-forming regions in the local universe fall into the "PS case," i.e., that superbubbles stall and fragment inside the disk and locally deposit almost all of their momentum. Higher-redshift observations suggest that z ∼ 2 SMGs exist in "PBO case" parameter space, i.e., superbubbles at z ∼ 2 are (always) able to drive outflows/fountains. The central regions of some local galaxies also appear to lie in the predicted "PBO case" region (e.g., NGC 4321, which has evidence of central star-formation-driven winds). The transition from high to low redshift galaxies, in terms of hosting pervasive outflows to only those in circumnuclear regions, appears driven by an evolution from high to low local gas fractions in star-forming regions according to this model.
• 2.  
  Observed intact H i bubble radii and velocities in local star-forming galaxies (from Bagetakos et al. 2011) are consistent with the ${\tilde{f}}_{\mathrm{gas}}$-Ω profile interpretations: most feedback driven bubbles in local galaxies should fragment inside the ISM, and that those able to break out originate from short dynamical time regions in the nuclear regions.
• 3.  
  A cluster formation model that includes a constant star formation efficiency per Toomre mass is effectively ruled out by the observational data (see Section 4.2), as this model predicts that high surface density regions (e.g., high-redshift star-forming clumps or low-redshift circumnuclear regions) would be unable to host superbubbles capable of breaking out and driving outflows/fountains.

In comparing to observations, we find that the clustering of SNe indeed has important implications for the local efficacy of star formation, and the evolution of galaxies more broadly across cosmic time. Future highly spatially resolved observations, capable of identifying and quantifying the properties of supernova-driven superbubbles, especially in dense molecular gas structures, should help to further constrain the effective strength of feedback under varying local galactic conditions and inform sub-grid models for feedback in cosmological galaxy simulations.