Any way the wind blows: quantifying superbubbles and their outflows in simulated galaxies across z ≈ 0-3

ABSTRACT

We present an investigation of clustered stellar feedback in the form of superbubbles identified within 11 galaxies from the FIRE-2 (Feedback in Realistic Environments) cosmological zoom-in simulation suite, at both cosmic noon (1 < z < 3) and in the local universe. We study the spatially resolved multiphase outflows that these supernovae drive, comparing our findings with recent theory and observations. These simulations consist of five Large Magellanic Cloud–mass galaxies and six Milky Way-mass progenitors (with a minimum baryonic particle mass of |$m_{\rm b.min} = 7100\,{\rm M}_{\odot }$|⁠). For all galaxies, we calculate the local and galaxy-averaged mass and energy-loading factors from the identified outflows. We also characterize the multiphase morphology and properties of the identified superbubbles, including the ‘shell’ of cool (⁠|$T\lt 10^5$| K) gas and break out of energetic hot (⁠|$T\gt 10^5$| K) gas when the shell bursts. We find that these simulations, regardless of redshift, have mass-loading factors and momentum fluxes in the cool gas that largely agree with recent observations. Lastly, we also investigate how methodological choices in measuring outflows can affect loading factors for galactic winds.

1 INTRODUCTION

Stellar feedback is one of the most important factors driving galaxy evolution. This feedback plays a critical role in the structure of the interstellar medium (ISM), as supernova (SN) explosions are a primary source of momentum, energy, and mass injection back into the ISM. Therefore, stellar feedback must be properly understood in order to accurately reproduce observed relationships and properties of galaxies, including the mass–metallicity relation (Tremonti et al. 2004; Davé, Finlator & Oppenheimer 2011; Ma et al. 2016; Wetzel et al. 2016; De Rossi et al. 2017; Kaviraj et al. 2017; Muratov et al. 2017; Chisholm, Tremonti & Leitherer 2018; Torrey et al. 2019; Porter et al. 2022; Bassini et al. 2024; Marszewski et al. 2024), turbulence (Faucher-Giguère, Quataert & Hopkins 2013; Kim & Basu 2013; Girichidis et al. 2016; Padoan et al. 2016; Orr et al. 2018, 2020; Burkhart 2021; Bieri et al. 2023), low star formation efficiencies (Hopkins et al. 2014; Burkhart 2018), stellar masses (Davé, Finlator & Oppenheimer 2012; Agertz & Kravtsov 2016; Wetzel et al. 2016), and the baryon cycle (Cen & Ostriker 2006; Muratov et al. 2015; Beckmann et al. 2017; Kim & Ostriker 2017; Anglés-Alcázar et al. 2017; Marinacci et al. 2018; Mitchell et al. 2020; Shin et al. 2023).

However, star formation and the resultant feedback do not occur uniformly within a galaxy. Instead, it occurs in clusters (Motte, Bontemps & Louvet 2018; Krumholz, McKee & Bland-Hawthorn 2019; Tacconi, Genzel & Sternberg 2020; Smith et al. 2021; Keller & Kruijssen 2022), and indeed Fielding, Quataert & Martizzi (2018) and Governato et al. (2010) note the necessity of such clustering in driving galactic winds. As a consequence of the proximity of formation of stars within a cluster, when core-collapse SNe from massive stars occur their shock fronts coalesce, and an even larger encompassing shock front known as a superbubble is born.

Superbubbles can be found in galaxies both at high redshift (z|$\ge$| 1) and in the local Universe (Taniguchi et al. 2001; Keller, Wadsley & Couchman 2015; Orr et al. 2022a; Watkins et al. 2023b). They typically consist of a hot, volume-filling component that sweeps up the surrounding ISM as it expands, creating a ‘hole’ in which very little cold gas resides. As the bubbles travel through the ISM, they capture cold molecular clouds that become entrained or homogenized within the hot and diffuse winds. This phenomenon results in the shell-like morphology of cold gas, referred to as the ‘cold cap’ of the superbubble in Orr et al. (2022b), around the energetic hot gas interior with evaporating cloudlet inclusions (Lancaster et al. 2021).

As a result of their importance, superbubbles have been widely studied for decades (Castor, McCray & Weaver 1975; Bregman 1980; Mac Low & McCray 1988; Ostriker & McKee 1988; Mac Low, McCray & Norman 1989; Koo & McKee 1992), but it is with recent computational methods and observational facilities that we truly become able to resolve the physics of superbubbles and their relationship with the baryon cycle and galaxy evolution. Notably, work enabled by JWST has begun to quantify the effects of superbubbles in nearby galaxies, such as NGC 628, in driving turbulence, sweeping up dense gas, and triggering additional star formation along the expanding shock fronts (Barnes et al. 2023; Mayya et al. 2023; Watkins et al. 2023a, b).

The disturbance of the ISM by superbubbles intimately connects galactic outflows and dense gas turbulence, as the successful breakout of bubbles from the galactic disc can drive some of the largest winds in galaxies (Muratov et al. 2015; Kim, Ostriker & Raileanu 2017; Fielding et al. 2018; Martizzi 2020). These multiphase winds, integral to the baryon cycle, are essential to the process of galaxy evolution (Strickland & Stevens 2000; Oppenheimer et al. 2010; Dalla Vecchia & Schaye 2012; Tumlinson, Peeples & Werk 2017; Faucher-Giguère & Oh 2023). Inflowing winds (meaning, reaccreting wind material) carry material necessary to sustain star formation and black hole growth, such as cold molecular gas, while new outflows from breakouts transport material into the circumgalactic medium (CGM). As a result, these winds (and the stellar feedback that drives them) are key to regulating star formation in galaxies (Oppenheimer et al. 2010; Ostriker, McKee & Leroy 2010; Hopkins et al. 2014; Hayward & Hopkins 2017; Guszejnov et al. 2022).

Outflow loading factors, often established in terms of mass (⁠|$\eta _{\rm M}$|⁠; Larson 1974; Veilleux, Cecil & Bland-Hawthorn 2005) and energy (⁠|$\eta _{\rm E}$|⁠; Larson 1974; Chevalier & Clegg 1985; Strickland & Heckman 2009), are used to quantify the properties of winds. These measures are evaluated as the respective outflowing quantity (e.g. mass or energy) normalized by the star formation rate, resulting in a dimensionless quantity that describes an efficiency relative to star formation. Loading factors have been frequently used in simulations to investigate and quantify galactic outflows and galaxy evolution (Strickland & Stevens 2000; Oppenheimer & Davé 2008; Dalla Vecchia & Schaye 2012; Muratov et al. 2015; Schaye et al. 2015; De Rossi et al. 2017; Nelson et al. 2019; Kim et al. 2020b; Mitchell et al. 2020; Pandya et al. 2021; Steinwandel et al. 2024), and while loading factors are more difficult to ascertain in observations, there still exists a sample to compare against theory (Martin 1999; Heckman et al. 2015; Chisholm et al. 2017; McQuinn, van Zee & Skillman 2019; Reichardt Chu et al. 2022; McPherson et al. 2023). However, measurements can vary across several orders of magnitude based on galaxy properties and dynamics. Several studies find that mass loading (⁠|$\eta _{\rm M}$|⁠) is dependent on current star formation surface density (⁠|$\Sigma _{\rm SFR}$|⁠) as |$\eta _{\rm M} \propto \Sigma _{\rm SFR}^{-0.44}$| (Li, Bryan & Ostriker 2017; Kim et al. 2020b; Li & Bryan 2020). Momentum-driven outflows, on the other hand, predict an even steeper dependence when outward velocity (⁠|${\it v}_{\rm out}$|⁠) is considered (⁠|${\it v}_{\rm out} \propto \Sigma _{\rm SFR}^2$|⁠; |${\it v}_{\rm out} \propto \rm SFR$|⁠; Murray, Ménard & Thompson 2011; Hopkins, Quataert & Murray 2012). Despite the variations in the dependence of outflows and loading factors on current and recent star formation, they have a well-studied relationship and aid in explaining the role of galactic winds in Large Magellanic Cloud (LMC)-mass galaxies (Fielding et al. 2018; McQuinn et al. 2019; Romano et al. 2023), and the mass–metallicity relation (Finlator & Davé 2008; Ma et al. 2016; Chisholm et al. 2018; Torrey et al. 2019; Bassini et al. 2024).

Simulations have become a powerful tool for studying superbubbles and the multiphase winds they can drive. However, the physical nature of these winds pose additional complications and restricts reliable studies to only simulations that can resolve a multiphase ISM. Another challenge for investigating superbubbles’ effect on their simulated host galaxies is that there is no standard method for identifying and quantifying galactic outflows. For example, some groups quantify outflows through fixed surfaces above and below galaxies (not unlike observational work) and take cuts on velocity relative to local escape velocity or the Bernoulli velocity (Kim et al. 2020b; Pandya et al. 2021), while others use particle tracking methods to differentiate outflows from fountains (Anglés-Alcázar et al. 2017; Hafen et al. 2019).

Orr et al. (2022a) developed an analytic model of clustered feedback from SNe, finding that the local gas fraction and dynamical time determine whether superbubbles broke out of the ISM (driving winds) or fragmented within the galaxy (driving turbulence). Orr et al. (2022b) compared this theory work with local observations from Barnes et al. (2023), Mayya et al. (2023), and Watkins et al. (2023b), finding general agreement. The work presented in this paper aims, in part, to compare these predictions to simulations in local galaxies and galaxies at cosmic noon.

In particular, this paper makes use of 11 galaxies from the FIRE-2 (Feedback in Realistic Environments¹; Hopkins et al. 2018) suite of cosmological zoom-in simulations, identifying and quantifying outflows from within these galaxies [five low-mass and six Milky Way (MW) mass] at both high redshift (z|$\sim$| 1–3) and the local universe (z|$\sim$| 0), analysing about 250 Myr of evolution at each redshift. These galaxies and cosmological epochs represent a diverse range over which we can study the launching of galactic winds and their properties, including superbubbles, as well as compare with recent observations.

We begin in Section 2 by detailing the FIRE-2 simulations and how we choose to define outflows within them. In Section 3, we present our results, including the wind properties and corresponding loading factors (Section 3.1), and their connection to superbubbles (Section 3.2). The physical implications of these results are further discussed in comparison to other simulations and observations in Section 4, and we briefly summarize our findings in Section 5.

2 SIMULATIONS AND ANALYSIS METHOD

We investigate superbubble feedback events in six MW/Andromeda-mass galaxies, and five Magellanic Cloud [Small Magellanic Cloud (SMC), Large Magellanic Cloud (LMC)]-mass galaxies (these masses all being achieved at roughly |$z\approx 0$|⁠), from the ‘standard physics’ FIRE-2 simulations introduced in Hopkins et al. (2018), which are publicly available (Wetzel et al. 2023). This work makes use of |$\sim$|10 snapshots of each galaxy near each integer redshift from zero to three (for a total of |$\sim$|40 snapshots per galaxy over its evolution), with an approximate time spacing between snapshots of |$\sim$|25 Myr, for a total elapsed time at each redshift of |$\sim$|250 Myr. Table 1 briefly summarizes the global gas and stellar mass properties of the simulations analysed here at each integer redshift.

The simulations analysed in this paper have baryonic particle masses of the order of |$m_{\rm b,min} = 7100$| M|$_\odot$|² and minimum adaptive force softening lengths <1 pc. Cooling in FIRE-2 is computed for gas temperatures |$T=10{-}10^{10}\,{\rm K}$|⁠. The suite of simulations includes a wide variety of heating and cooling physics, including free–free, photoionization/recombination, Compton, photoelectric, metal-line, molecular, fine-structure, and dust collisional processes. In particular, metal-line cooling is noted by Hopkins et al. (2018) to be particularly important for superbubbles. The gas softening lengths are adaptive and we note that the effective radius of gas elements at the minimum density of star formation is about 7 pc (⁠|$n=1000\,\rm cm^{-3}$|⁠; see table 3 of Hopkins et al. 2018), with the densest structures having shorter smoothing lengths down to the sub-pc minimum.

In the FIRE-2 simulations, stars form on a free-fall time in gas that is dense (⁠|$n \gt 10^3$| cm|$^{-3}$|⁠), molecular (according to the Krumholz, Leroy & McKee 2011 methodology), self-gravitating (viral parameter |$\alpha _{\rm vir} \lt 1$|⁠), and Jeans-unstable. Star particles are considered single stellar populations with defined age, metallicity, and mass.

The FIRE-2 ‘standard physics’ incorporates feedback mechanisms from SNe, stellar winds from OB-type stars/asymptotic giant branch (AGB) stars, photoionization and photoelectric heating, and radiation pressure. This suite excludes active galactic nucleus (AGN), cosmic rays, and additional magnetohydrodynamic physics, although other studies within the broader FIRE-2 project have explored these ‘extended physics’ elements (Anglés-Alcázar et al. 2017; Chan et al. 2019; Su et al. 2019). For comprehensive details on the ‘standard’ physics and their application, see Hopkins et al. (2018). Of particular importance to this current study of superbubbles/SN-driven winds are the core-collapse and Type Ia SN rates, derived from starburst99 (Leitherer et al. 1999) and Mannucci, Della Valle & Panagia (2006), respectively.

We generate maps of the gas, stellar, and star formation rate properties from the snapshots using the same methods as Orr et al. (2018, 2020), projecting the galaxies face-on (or edge-on) using the angular momentum of the star particles within the stellar half-mass radius, and binning star particles and gas cells into square pixels with side-lengths (i.e. ‘pixel sizes’ |$l_{\rm pix}$|⁠) 750 pc. The maps are 30 kpc on a side, and integrate gas and stars within |$\pm 15$| kpc of the galactic mid-plane. Cold gas structures are well contained within a single one of these pixels, and the most diffuse hot gas cells are marginally resolved by these pixels down to densities of |$\sim 10^{-3}$| cm|$^{-3}$|⁠. The length-scale is comparable with the resolvable scale of JWST observations at |$z\sim 1$| (Boker & Arribas 2022).

We generate a proxy for observational measures of recent SFRs by calculating the 40-Myr-averaged SFR. We do this by summing the mass of star particles with ages less than 40 Myr, and correcting for mass-loss from stellar winds and evolutionary effects using predictions from starburst99 (Leitherer et al. 1999) to serve as correction factors. FIRE-2 also uses these predictions to tabulate all feedback quantities, determine the ionizing photon budget for all star particles, and provides estimates for approximate mass-loss rates. For further information on the role of starburst99 in FIRE-2, we direct the reader to section 2.5 of Hopkins et al. (2018). The 40 Myr time interval was chosen for its approximate correspondence with the time-scales traced by continuum ultraviolet (UV)-derived SFRs (Lee et al. 2009).

To calculate the outflow mass, momentum, and energy fluxes in a snapshot, we project the gas cells on to planes a fixed height above/below the galaxy and then calculate the flux quantities through that surface in a Cartesian grid of 750 pc pixels using the following definitions:

and

where i is summing over all the gas cells whose kernel overlaps with the pixel surface, |$\rho$| is the gas element density evaluated at the pixel surface, |${\it v}_z$| is the gas velocity normal to the pixel surface, |$v^2$| is the square of all velocity components,³|$\gamma =5/3$| is the adiabatic index (i.e. heat capacity ratio) for a monatomic gas, and u is the specific internal energy of the gas element. For high-redshift snapshots (⁠|$z = 1{-}3$|⁠), we select a height above/below the main galaxy body of 0.05|$R_{\rm vir}$|⁠, and for the low-redshift snapshots (⁠|$z \approx 0$|⁠) we select a height above/below the galaxy of twice the gas scale height |$2H$|⁠. We note that while this is a different height than 0.05|$R_{\rm vir}$|⁠, our choice of 2H at late times maintains the qualitative spirit of investigating the height at which winds are launched from the ISM, where the direct connection between young stellar populations in the ISM and wind loadings is most measurable. We note that where outflows are defined can have an effect on the resultant loading factors, and we briefly analyse how different choices for the height of the flux surface affect our results in Appendix  A.

2.1 Characterizing the individual simulations

We find it important to briefly review the galaxies’ basic properties and behaviour. While we use 11 FIRE-2 galaxies in this study, the majority of examples of specific superbubble events shown in the Figures will be from two galaxies for conciseness: one from the low-mass m11 galaxies, m11d, and one from the MW-mass galaxies, m12b. This allows us to showcase superbubbles and wind properties of both galaxy types, in addition to bubbles occurring at different redshifts. Table 1 presents the stellar masses, gas masses, and gas fractions (⁠|$f_{\rm gas} = M_{\rm gas}/(M_{\rm gas}+M_{\rm \star })$| at the radius within which 90 per cent of the galaxy snapshot’s stellar mass is contained) of the simulated galaxies analysed in this paper in each of the redshift bins from |$z = 3{-}0$|⁠.

Our sample of low-mass galaxies (m11s) shows remarkable diversity in morphologies and modes of star formation for five galaxies of similar mass. Only one LMC–mass galaxy, m11h, forms a disc, while the others are irregular galaxies with small starburst events. They generally remain high in gas fraction (⁠|$f_{\rm gas} \gtrsim 0.5$|⁠) for their entire evolution down to |$z =0$|⁠. The top panels of Fig. 1 illustrate the spatial distribution of all gas throughout galaxy m11d, with pixel size 750pc: This galaxy shows little distinct morphological features, though signs of strong stellar feedback are present at |$z\sim 3$|⁠.

The larger galaxies in our sample (the m12s) all form discs by |$z=0$| and have lower gas fractions than the m11s even by |$z =3$|⁠. The m12s consume much more of their gas in star formation, and their gas fractions fall earlier in cosmic time, and more dramatically, than the m11s resulting in gas fractions of |$f_{\rm gas} \approx 0.1{-}0.3$|⁠. The |$z=0$| panels of the bottom section of Fig. 1 shows the discy morphology of m12b at |$z=0$|⁠: we can clearly see spiral arms |$z\sim 0$|⁠, and the edge-on view shows a clear disc that was not evident at higher redshift.

One of the most distinct features present in all 11 galaxies, across all redshifts, is the frequent stellar feedback. It is well established that many FIRE galaxies are relatively bursty in star formation before the formation of discs (Sparre et al. 2017; Faucher-Giguère 2018; Orr et al. 2018, 2020, 2021; Stern et al. 2021; Gurvich et al. 2023; Hopkins et al. 2023; Sun et al. 2023). Specifically, El-Badry et al. (2017), Anglés-Alcázar et al. (2017), and Sparre et al. (2017) note that the resolved ISM and stellar feedback physics in FIRE-2 gives rise to ‘breathing modes’ of star formation that often continue up until disc formation (z|$\approx 0.4{-}0.7$| for MW-mass galaxies and m11h).

3 RESULTS

3.1 Mass and energy loadings

As previously mentioned, loading factors are useful for the study of galactic outflows due to their ability to relate outgoing mass and energy with star formation. In line with convention, we calculate the mass-loading factor |$\eta _{\rm M}$| as

and the energy-loading factor |$\eta _{\rm E}$| as

where (⁠|$E_{\rm SN}$|/100 M|$_\odot$|⁠) represents the mechanical energy injection rate per mass of stars formed, which is |$10^{51}$| erg per 100 M|$_{\odot }$|⁠. Both loading factors are dimensionless measures of outflows, normalized by the contribution from star formation.

Here, we use the 40 Myr-averaged SFR to calculate |$\eta _{\rm M}$| and |$\eta _{\rm E}$|⁠, as this SFR can be inferred from UV observations, and represents time-scales long enough to trace injection of energy from core-collapse SNe from young clusters. We calculate values for the loading factors based on the overall outflows (gas with |${\it v}_{\rm out} \gt 0$|⁠). We also compare the effects of imposing a velocity cut (three times the neutral gas velocity dispersion in the disc: |${\it v}_{\rm cut} \gt 3\sigma _{\rm neut,z}$|⁠) on our local loading factors, which represent more conservatively estimated outflows, in Fig. 6.

In presenting our loading factors, we compute these measures both locally in the galaxies (per 750 pc pixels) and globally (averaged in each galaxy snapshot). Global values are divided by the area of the galaxy, assuming a circular shape (⁠|$\pi r^2$|⁠) approximated using the radius at which half of the SFR resides (cf. a galactic H|$\alpha$| radius). Local (pixel) values are divided by |$\ell ^2$|⁠, where |$\ell$| is the pixel size of 750 pc.

3.1.1 Local mass loading

Fig. 2 shows the spatially resolved distribution of |$\eta _{\rm M}$| with star formation rate surface density, with the top row showing the cold–warm (⁠|$T\lt 10^5$| K) gas and the bottom row the hot (⁠|$T\gt 10^5$| K) gas phase. Each column represents the ten snapshots in each redshift range analysed here (⁠|$z\sim 3{-}0$|⁠). Points represent the median |$\eta _{\rm M}$| value within 0.75 dex wide bins of SFR surface density (where each bin must have at least 50 pixels), with corresponding 1|$\sigma$| error bars, and lines are the best-fit of these points. We also demarcate the region in which thermal CGM motions may dominate our measured ‘outflow’ fluxes, see Section 3.1.5 for details.

We show this distribution for both phases of the gas (hot and cold–warm) across the redshifts, and also separate our galaxy populations into the m11s (blue) and m12s (orange). For each galaxy sample, we plot a best fit of the median |$\eta _{\rm M}$| values in each redshift panel. Equations for these lines of best fit can be found in Table 2. We also plot predictions from the simulations of Kim et al. (2020b) and Steinwandel et al. (2024).

Comparing the distributions of SMC–mass galaxies (m11s) versus MW-mass progenitors (m12s), we can see that at all times the two have very similar distributions of |$\eta _{\rm M}$|⁠, with the m12s hosting higher SFRs at cosmic noon, resulting in the fits extending out to |$\dot{\Sigma }_{\star }=10^1$| M|$_\odot$| pc|$^{-2}$|⁠, while the m11s only reach about |$10^{-0.5}$| M|$_\odot$| pc|$^{-2}$|⁠. Both galaxy samples maintain identical slopes of |$\eta _{\rm M}$|⁠. The cold–warm gas shows little signs of evolution with redshift, while the hot gas mass loadings show a subtle decrease at z|$\sim$|0 compared to z|$\sim$|1–3 (this is muddled somewhat by the change in the altitude at which the loading is measured).

Our spatially resolved cold–warm gas mass-loadings are between the predictions from the highly resolved isolated galaxy simulations of Steinwandel et al. (2024), and the tall-box TIGRESS simulations from Kim, Choi & Kim (2020a). It appears that the cooler gas in all of our FIRE-2 galaxies, at all redshifts, is carrying slightly more outflowing mass than the hot gas phase, consistent with other FIRE-2 measurements by Pandya et al. (2021).

The spatially resolved hot gas |$\eta _{\rm M}$| is very similar to the cold–warm gas in form with a similar power-law slope, though with an overall slightly lower normalization. This is in contrast with predictions of simulations by Kim et al. (2020a) and Steinwandel et al. (2024), who both predict a nearly flat relationship in the hot gas between |$\eta _{\rm M}$| and |$\dot{\Sigma }_{\star }$|⁠. We, however, see little differences between the cold–warm gas and hot gas.

3.1.2 Global mass loading

We also calculate |$\eta _{\rm M}$| globally, summing the outflows and star formation over the entire galaxy at each point in time to calculate the loading factors. Fig. 3 presents the galaxy averaged mass loadings as Fig. 2, and include observations from McQuinn et al. (2019) plotted as x’s in the upper right panel, and simulation snapshot points with that include an identified superbubble outlined in black. Fit lines for global values of the loading factors can be found in Table 3.

With the global mass-loading values, we see an opposite slope as compared to the local calculations in Fig. 2, and the two galaxy groups now occupy different spaces in |$\dot{\Sigma }_{\star }$| and |$\eta _{\rm M}$|⁠. The SMC-mass m11 galaxies have lower star formation rate surface densities, as might be expected. And they also appear to have higher mass-loading values at all redshifts relative to their global SFRs, speaking to the ability of winds escaping these lower potential galaxies. We find good agreement with observations from McQuinn et al. (2019), though the FIRE-2 m11 galaxies have higher SFRs as a sample. As with the local scale, on the galaxy-scale mass-loading is dominated by the cold–warm gas. The cool gas mass loading is |$\sim$|dex higher than the hot gas, across redshift.

When it comes to the presence of superbubbles in these galaxy-averaged loading factors, the m11 z|$\sim$| 0 superbubble snapshots tend to lie at the higher end of |$\dot{\Sigma }_{\star }$|⁠, though with lower |$\eta _{\rm M}$| than the best fits (see second column of Fig. 3). The superbubbles in the m12s show higher mass-loadings in the hot gas at high redshift than the best-fitting line (points with grey outlines in lower left panel of Fig. 3), but they appear to be more scattered in the cooler gas (upper left panel). We note that no significant superbubbles were identified in the m12 |$z\approx 0$| snapshots, and the bubble snapshots we have identified also include one snapshot before and after the bubble itself is identified, which may affect results. In particular, the m11 superbubble atz|$\sim$| 0 shows that each snapshot of the bubble (before, during, and after) has significantly different loading factors in the hot gas (see lower right panels of Figs 3 and 5).

3.1.3 Local energy loading

Fig. 4 displays the same information as Fig. 2, with the exception that we now plot the local energy-loading factor, |$\eta _{\rm E}$|⁠, instead of |$\eta _{\rm M}$|⁠. As with |$\eta _{\rm M}$|⁠, the fit values for local |$\eta _{\rm E}$| can be found in Table 2.

The energy-loading factors of the low-mass galaxies and MW-mass progenitors appear to differ more than the mass loading factors. The MW-mass m12s, at all redshifts, have energy-loading factors that are nearly always larger than the m11s, reaching up to an order of magnitude difference in the cooler gas at z|$\approx$|0 (see upper right panel of Fig. 4). This difference is not quite as large at higher redshifts. While it may be expected that low-mass galaxies (e.g. m11s) have higher energy-loading factors due to their shallower potentials/more ‘efficient’ outflows, we note that measuring outflows in the same way for both galaxy types can result in an inherent bias: larger galaxies like the m12s may need to have significantly more energy in their outflows in order to expel mass to the same relative height (⁠|$\rm \pm 0.05 {\it R}_{vir}$| or 2H) due to their larger gravitational potentials.

When comparing to Kim et al. (2020a) and Steinwandel et al. (2024), we find that our energy loadings in the cold–warm gas generally are higher, though there is significantly less tension when considering the m11s. We find opposite slopes in the hot gas as compared to their simulations: while Kim et al. (2020a) and Steinwandel et al. (2024) show positive slopes of |$\dot{\Sigma }_{\star }$| and |$\eta _{\rm E}$|⁠, we maintain negative slopes, as in the cold–warm gas. One possible explanation for this is that when calculating loading factors, we use the 40 Myr-averaged star formation rates. Some measurements of energy-loading may be biased due to this: any amount of detected ‘outflowing gas’ that was launched by the tail end of a starburst |$\sim$|40 Myr ago could dramatically inflate loading factor values, whereas there is not a clear way in which energy loadings could be diluted to an extreme. This caveat may affect the resulting |$\eta _{\rm E}$| or |$\eta _{\rm M}$| distributions.

3.1.4 Global energy loading

As with |$\eta _{\rm M}$|⁠, we also show the global (galaxy averaged) distribution of the energy-loading factors for this sample in Fig. 5, with fits available in Table 5.

Distributions in the global |$\eta _{\rm E}$| are very similar to those from Fig. 3. The MW-mass progenitors appear to have higher global SFR surface densities than the low-mass m11s, and energy-loading factors at high redshift are similar (albeit shifted in SFR). At high redshift, the energy loading is dominated by the hot gas, which matches theoretical predictions. At low redshift, however, the cooler gas appears to have slightly more energy, especially in the m11s. This may well be related, however, to the lack of significant large-scale outflows from the m12s, and generally smooth (in time) SFRs.

3.1.5 Contributions from thermal motions and/or CGM gas

To conservatively estimate the contribution of thermal motions in the inner CGM/halo of the simulations to our mass and energy flux measurements (essentially a resolution limit), we calculate here estimates for the mass and energy flux of a single gas element with representative temperatures and densities moving with thermal velocity across the flux measurement surface, yielding essentially minimum |$\eta _{\rm M}$| and |$\eta _{\rm E}$| values as a function of |$\dot{\Sigma }_{\star }$| for the cold–warm and hot gases. To do so, we calculate fiducial values using the minimum resolvable measurements of these particular FIRE-2 simulations.

Estimating the mass and energy fluxes, for a single particle at the minimum baryonic particle mass, given its effective size (calculated with |$m\approx \rho h^3$|⁠) we have: |$\dot{M} = h^2 \rho v$| and |$\dot{E} = h^2 \rho v (v^2 + c_{\rm s}^2/(\gamma -1))$|⁠, where we have taken the standard sound speed definition |$c_{\rm s}^2 = (\gamma -1) u$| and used equations (1) and 3. Assuming that the velocity of the particle is the thermal sound speed (see Fig. 6 for a distribution of local outflow velocities compared to the sound speed, these fluxes simplify to |$\dot{M} = h^2 \rho c_{\rm s}$| and |$\dot{E} = h^2 \rho c_{\rm s}^3 \gamma /(\gamma -1)$|⁠. For fiducial values, we take the minimum baryonic particle mass in these FIRE-2 simulations, |$m_{\rm b, min}=7100\,{\rm M}_{\odot}$|⁠, and the velocity to be the sound speed of the gas (calculated using the standard definition above, resulting in for the cool gas, |$c_{\rm s} \approx 12$| km s^-1 assuming |$\mu =0.59$| and |$T=10^4$| K, while the hot gas has |$c_{\rm s} \approx 35$| km s^-1 around |$T=10^5$| K). Then, the smoothing length h is found by assuming densities of |$10^{-2}$| and |$10^{-3}$| cm|$^{-3}$| for the cool and hot gases, respectively, resulting in |$h_{\rm cool}\approx 300$| pc and |$h_{\rm hot}\approx 650$| pc.

The estimated ‘thermal’ loading factors then scale by star formation surface density as

and

with both the numerator and |$\dot{\Sigma }_{\star }$| in units of M|$_{\odot }$| yr|$^{-1}$| kpc|$^{-2}$|⁠.

These ‘thermal’ loading factors are visible as the shaded regions in Figs 2–5. The grey-shaded region essentially representing where our outflow measurements may be dominated by thermal motions of a single gas element at the simulation resolution scale in the inner CGM, as these particles are travelling at or somewhat below the sound speed.

Comparing this thermal limit with our spatially resolved distributions, we see that as much as half of the hot |$\eta _{\rm M}$| values at |$z\approx 0$| (see lower right panel of Fig. 2) fall in this potentially unresolved thermal regime. The cool gas mass loadings are generally above this cutoff (i.e. it is highly supersonic), but at |$z\approx 0$| the low SFR end of the cold–warm gas |$\eta _{\rm M}$| distribution abuts the thermal regime.

Concerning the spatially resolved energy-loading factors |$\eta _{\rm E}$|⁠, the pixels are generally well above the thermal regime, though again the hot gas at |$z\approx 0$| has some significant overlap. This highlights the importance of resolving this diffuse phase of galactic winds/the CGM, especially in simulations with weak winds (like these FIRE m12 runs at |$z\approx 0$|⁠). We show in Fig. 6 that our velocity cut of |$3\sigma _{z}$| does not significantly help to separate our loading distributions any further from this ‘thermal’ regime. This suggests that a (larger) outflow velocity cut does not allow us to any more cleanly interpret our results, rather that we appear to be at least marginally resolving the outflow dynamics in these simulations.

3.2 Superbubbles within the simulations

3.2.1 Identifying superbubbles

We visually identified superbubbles within the simulations using the total gas distributions (see Fig. 1) from all snapshots in our sample (oriented face-on). If a clear bubble-like morphology was present in the gas, we then examined the snapshots directly before and after the ‘bubble’ looking for additional tracers of superbubbles, similar to observational methods in Watkins et al. (2023b). These include the ‘shell’ of cooler gas, as well as a diffuse hot gas component within the shell. In addition to these quantities, we used the 40 Myr-averaged SFR to trace whether an ‘active’ star cluster was at the centre of this cooler gas shell (producing clustered SNe to drive the superbubble). If all of these conditions were met, we identified the event as a likely superbubble.

We conservatively identify a total of 9 primary superbubbles occurring throughout these 11 galaxies and four redshift ranges, with four of the superbubbles belonging to the m11 galaxies and five to the m12s. Certainly many more superbubbles occur between snapshots, or are not as clean-cut, so this is a lower limit to their occurrence rate in these simulations, and we highlight these as archetypal examples of this feedback mode. For the sake of conciseness, we present results in the main text for two specific superbubbles (the remainder can be found in Appendix  B): one in m11d at |$z\sim 0$|⁠, and one in m12b at |$z\sim 1$|⁠. These two particular superbubbles allow us to showcase results in one low-mass galaxy and one MW-mass progenitor, as well as at both low and high redshift.

Fig. 7 shows the cold–warm and hot gas surface densities, and the star formation rate surface density, for the two superbubble events of interest (left panel: m11d, right panel: m12b). Both galaxies are oriented ‘face-on’, however, lacking a disc this distinction is not terribly important due to the isotropic gas distribution in high-redshift and low-mass galaxies. We note that viewing angle may have an impact on galaxies when a disc is present, discussed further in Section 4.

On the left side, m11d’s superbubble is shown, with the first column being the snapshot before the superbubble (about 22 Myr beforehand), the middle column showing the first snapshot of the shell of cooler gas expanding, and the third column as the snapshot immediately afterwards. The three columns together show a total span of about 44 Myr in time from the first column’s snapshot. In the top row, we have the surface density of the cold–warm (⁠|$T\lt 10^5$| K) gas, where the bubble itself is most visible. The middle row shows the 40 Myr-averaged star formation rate surface density, and the bottom row is the diffuse hot (⁠|$T\gt 10^5$| K) gas. The right side of Fig. 7 shows the same information, but for galaxy m12b, spanning a total of about 50 Myr in time. These panels of m12b depict a particularly stunning example of what one might expect from a superbubble: stellar feedback from the central star cluster couples together, creating the signature bubble shape in the cold–warm gas, which then breaks open and fragments, releasing the hot gas component. The bubble in m11d appears to reach about 6 kpc in diameter before it fragments in the third snapshot, while the bubble in m12b is much smaller at 3–4 kpc in diameter.

3.2.2 Superbubble winds: mass and energy fluxes and the connection to star formation

Fig. 8 depicts the total outflowing mass and energy flux from each galaxy in Fig. 7, m11d and m12b, during the redshift range for which the bubble was identified. We plot both the hot and cold–warm gas, and identify the time frame of the snapshots shown in Fig. 7 with a shaded region |$\pm 1$| snapshot. Lighter shaded regions, if present, show smaller or secondary instances of stellar feedback.

In both m11d and m12b, we see a sharp rise in the outflowing mass in all gas phases during the identified superbubble events. In m11d at redshift zero, during the snapshots of the superbubble, the total hot gas mass flux increases by more than an order of magnitude, while the cold–warm gas increases about 0.5 dex. The outflows steadily decrease after the superbubble’s fragmentation (right edge of shaded region). m12b shows a similar story: The hot and cold–warm mass fluxes increase by 0.5–1.5 dex, then (on average) decreasing, though we do see another increase in the fluxes during this time span, suggesting a secondary event of strong stellar feedback. For both m11d and m12b, we note that the cooler phase of the gas carries a majority of the mass relative to the hot gas.

As with the mass flux, we see an increase in the energy flux during the identified superbubble events, followed by a significant decrease. In m11d, the phase that seems to experience the most dramatic changes is the hot gas, rising up three orders of magnitude before decreasing, while in m12b the energy flux of both phases is almost equal. In m12b, we can see that the hot gas energy flux decreases by about two orders of magnitude some time after bubble fragmentation. In both events shown here, it is apparent that the outflowing energy flux of hot gas reaches a peak during the identified superbubble.

In both galaxies, we see a significant increase in the star formation rate as the superbubble occurs, as much as an order of magnitude for m11d, followed almost immediately by a decrease after the bubble’s fragmentation and release of hot and diffuse gas and energy. m11d shows a much clearer relationship, as the galaxy as a whole is experiencing less disruption at low z than m12b near |$z=1$|⁠. Tellingly, for every superbubble identified within the galaxies of this sample (see Appendix  B), we observe a direct relationship with the star formation rate (dashed black line) and the mass/energy fluxes. Given the 40 Myr SFR averaging window, and the steep rise in SFRs with the bubble expansion, it is clear that the majority of star formation associated with the superbubble events occurs in |$\lesssim 25$| Myr as the bubble is first expanding. This supports recent findings from JWST of a temporary increase in star formation during superbubbles in NGC 628 (Barnes et al. 2023; Watkins et al. 2023b), but it is not entirely clear that ‘triggered star formation’ is a dominant scenario associated with superbubbles as opposed to a single large starburst.

3.3 What are the conditions for superbubble break out?

Orr et al. (2022b) presented a predictive model for whether a superbubble will ‘break out’ of the galactic disc (powered breakout, PBO), or stall within (powered stall, PS), based on parameters of the ISM, namely the local gas fraction |$f_{\rm gas}$| and inverse dynamical time |$\rm \Omega ={\it v}_{\rm c}/R$|⁠. They then compare this theory with both simulations and observations. Orr et al. (2022b) make several predictions regarding the balance of momentum injection into the ISM versus CGM as a function of ISM properties, arguing that perhaps as much as |$\sim 60~{{\ \rm per\ cent}}$| of SN momentum will go into the outflow, as opposed to the ISM (with observational evidence corroborating this by Reichardt Chu et al. 2022). None of the model assumptions require a disc environment, and can therefore be tested by the range of galaxy morphologies/masses/sizes we use in this study (see section 7.4 of Orr et al. 2022a).

Fig. 9 shows our results for spatially resolved local (⁠|$z\approx 0$|⁠) momentum fluxes across our sample, plotted in the |$\Omega$|–|$f_{\rm gas}$| phase space used by Orr et al. (2022a, b). We do not distinguish between our SMC/LMC-mass galaxies (m11s) and MW-mass progenitors (m12s) here, instead organizing our results into bins of stellar mass. The shaded region shows the distribution of spatially resolved (750-pc scale) outflows at |$z\approx 0$|⁠, coloured by the average momentum flux value. A solid black line denotes the boundary line between powered superbubble breakout (PBO; right side of the boundary line) and PS (left side of the boundary line) within the disc (see equation 9 of Orr et al. 2022a). For comparison, spatially resolved observations from a starbursting disc galaxy presented in Reichardt Chu et al. (2022) are plotted as smaller points with black backgrounds in the upper right panel of Fig. 9 (cold–warm gas and|$M_\star \gt 10^{10}$| M|$_{\odot }$|⁠).

The evolution in gas fractions and inverse dynamical times with stellar mass is evident, necessitating the separation of stellar mass into various bins. As stellar mass increases, |$\Omega$| increases, and |$f_{\rm gas}$| decreases. Similar momentum fluxes are seen in both the hot and cold–warm gas for stellar masses |$\lt 10^{10}$| M|$_\odot$|⁠. However, in the galaxies with larger stellar masses, the hot gas carries |$\gtrsim 0.5$| dex more momentum than the cold–warm gas. This may be explained by winds in non-disc (and lower mass) environments being more mixed, and thus carrying similar momenta. However, for the |$\gt 10^{10}$| M|$_\odot$| galaxies, the hot gas carrying more momentum than the cooler gas is suggestive of more ordered multiphase winds that are less homogeneous.

Comparing our spatially resolved cold–warm gas momentum fluxes with observations by Reichardt Chu et al. (2022), their values coloured on the same scale in the Figure, in the range of about |$10^{-1}{-}10^{0.50}$||$\rm\, {\rm M}_{\odot} \,km s^{-1}\,yr^{-1}$|⁠, it is clear that their observed galaxy IRAS08 is significantly more gas rich than the FIRE-2 m12s in its centre and hosts a far stronger outflow. However, there is agreement where our distributions overlap and generally a similar trend. More extreme (in gas fraction, and presumably SFR) FIRE-2 discs would make for an excellent test of both the FIRE-2 model and the predictions of Orr et al. (2022b).

In Fig. 10, we show global (galaxy averaged) outflow momentum fluxes in |$f_{\rm gas}$|–|$\Omega$| space from |$z\sim 1{-}3$| snapshots. As with Figs 3 and 5, we again mark the superbubble snapshots with black outlines. Similar trends in gas fraction and dynamical time are seen, with decreasing |${\it f}_{\rm gas}$| and increasing |$\Omega$| with stellar mass and redshift. Momentum flux values clearly also increase with increasing stellar mass than with redshift alone.

Most points lie in the PBO region from Orr et al. (2022b), though we note that not every superbubble snapshot does. This may occur due to the ‘superbubble snapshots’ including a snapshot before and after the bubble is first identified, and gas could be highly depleted following a superbubble disrupting the ISM. These superbubbles only seem to have among the highest momentum fluxes in each respective mass bin at the lowest stellar masses of |$M_\star \lt 10^{9} {\rm \,\, M_\odot }$|⁠. At the highest stellar mass bins (⁠|$M_\star \gt 10^{9.5} {\rm \,\, M_\odot }$|⁠), they have among the lowest gas fractions.

4 DISCUSSION

It has been well established that stellar feedback in the form of superbubbles drives some of the largest galactic-scale winds. However, the exact detailings and predominant physics involved in this process have been a topic of debate, particularly relating to the balance of regulating star formation in the ISM versus driving winds, the phase structure of the CGM, and various outflow properties.

Both Kim et al. (2017) and Fielding et al. (2018) highlighted the importance of spatial and temporal clustering of SNe in driving superbubble winds, without which the SN remnants would simply be expected to merge with the surrounding turbulent ISM. The clustering, tied to the mass of the young star cluster, led to direct predictions by Fielding et al. (2018) of loading factors for these winds. They predicted that energetic winds should have |$\eta _{\rm M} \approx 0.5{-}1$| at |$z\sim 2$|⁠, where the mass ejected by stellar feedback is of the order of the star formation rate, reminiscent of a ‘bathtub model’ (Dekel & Mandelker 2014; Belfiore et al. 2019). These predictions are similar to what we see in our FIRE-2 galaxies at cosmic noon (evident in Figs 2 and 3), particularly in the cooler gas.

Furthermore, Fielding et al. (2018) suggested that high mass-loading factors (⁠|$\eta _{\rm M} \gg 1$|⁠) in low-mass galaxies, similar to what we find at high redshift for the m11s, are often halo-scale quantities, as efficient venting of superbubble winds leads to less accretion on to the galaxy from the CGM, increasing the value of |$\eta _{\rm M}$| when averaged over scales comparable to the halo virial radius. Several studies have also suggested that hot outflows transport a majority of the energy from stellar feedback into the CGM (Chevalier & Clegg 1985; Li et al. 2017; Kim et al. 2020b; Li & Bryan 2020; Fielding & Bryan 2022). The argument follows that these hot energetic outflows deposit their heat in CGM, subsequently regulating star formation in the galaxy by preventing cold gas accretion (Ostriker et al. 2010; Hopkins et al. 2014; Hayward & Hopkins 2017; Li & Bryan 2020; Fielding & Bryan 2022). By this logic, we might expect that galaxies with more efficient hot and energetic winds, such as low-mass galaxies with their shallower potential wells, will experience stronger star formation rate quenching following a large superbubble feedback event (Muratov et al. 2015; Pandya et al. 2021).

We find this result in our simulations. As mentioned in Section 3.2.2, a clear correlation exists between extremely energetic hot winds from a galaxy, usually occurring during a superbubble we have identified, and the 40 Myr-averaged star formation rate. Appendix  B contains a more complete picture of this result, showing our results for all the superbubbles, we identify in this study. We note that the simulated galaxies, regardless of mass, experience a drop in star formation of about an order of magnitude following the superbubble’s breakout and shell fragmentation, relative to their pre-superbubble SFRs. The drop in SFRs does appear to be more pronounced in the SMC-mass galaxies (see e.g. m11d and m11e in Fig. B1), regardless of redshift. And so, when the gravitational potential is lower winds appear to be far more effective at escaping the galaxy and suppressing star formation in FIRE-2 galaxies.

Relatedly, Orr et al. (2022a), using observational data, generally substantiated the suggestion by Orr et al. (2022b) that superbubble feedback at higher redshifts (⁠|$z\ge 1$|⁠) regulate the ISM differently: they predicted that superbubbles at |$z\sim 2$| should always drive true outflows and fountains, injecting most of the SN energy/momentum into the CGM. Our results here show a clear evolution in gas fraction and orbital time over redshift – this is largely linked to the corresponding evolution in stellar mass (see Fig. 9) – overall finding that gas fraction increases with increasing redshift, while the orbital time decreases. Nearly all of our galaxy snapshots, at all redshifts, lie within the ‘powered breakout’ region from Orr et al. (2022b), indicating that these superbubbles should indeed drive galactic fountains and outflows if they occur. Though when we spatially resolve the inner regions of the m12s (⁠|$M_\star \gt 10^{10}$| M|$_\odot$|⁠), Orr et al. (2022b) would predict that those regions should not host significant, powered superbubble breakouts. To wit, this aligns with our findings, whereas for generally all of the spatially resolved regions of all of the lower mass galaxies, Orr et al. (2022b) predict superbubble breakout, i.e. none of the lower mass galaxies should be expected to contain their SN feedback, just as we observe in our study.

4.1 Comparison to other simulations

Throughout the paper we have made comparisons with results from both the tall-box TIGRESS simulations (Kim et al. 2020b), and those of Steinwandel et al. (2024). Their works both study SN-driven outflows in extremely high resolution, explicitly modelling a self-consistent ISM, feedback, and star formation at small scales: 2–8 pc in spatial resolution in a 500 pc tall-box disc patch in TIGRESS, and at |$\sim$|4 M|$_{\odot }$|⁠, |$\sim$|1 pc resolution in an isolated LMC-mass galaxy in Steinwandel et al. (2024). However, we caution the reader against too close a comparison with these simulations as the simulation methods are vastly different, with neither of their simulations including a realistic CGM into which the outflows are launched. In addition, the number of caveats in calculations can make comparison difficult, as no ‘gold standard’ method exists for calculating outflows. We use slightly different temperature cuts on the gas, spatial resolution (750 pc), a 7100 M|$_{\odot }$| mass resolution, and use a more diverse sample of galaxy types and redshifts.

Our spatially resolved |$\eta _{\rm M}$| values for the cold–warm gas in Fig. 2 fall between the of mass-loading values predicted by Kim et al. (2020b) and Steinwandel et al. (2024). Interestingly, the FIRE-2 galaxies we analysed have near identical dependencies for both gas phases on SFR, though with a lower normalization for the hot gas (it does not carry much of the mass). And although Steinwandel et al. (2024) and Kim et al. (2020b) use similar definitions for hot gas, we find a steeper dependence on SFR for mass loading of the hot phase, which may be related to our ability to resolve the hot phase on |$\sim$|kpc scale (motion at roughly the sound speed could dominate our hot gas fluxes). Indeed, from Fig. 6, we see that much of the gas is supersonic, and outflow velocities tend to decrease with increasing cosmic time. Furthermore, the hot gas typically travels at larger velocities than the cool–warm phase, consistent with previous outflow predictions.

When considering the spatially resolved energy-loading factors, we find considerably more tension between our results and those of Kim et al. (2020b) and Steinwandel et al. (2024). Generally speaking we find significantly more energy in the cold–warm gas phase, except for winds launched from the most star forming of regions (top panels in Fig. 4), though the SMC-mass galaxies (m11s; blue in Fig. 4) are closer in magnitude. In contrast with these previous simulation works, we see little differences in the slopes of |$\dot{\Sigma }_{\star }$| versus |$\eta _{\rm E}$| between the cold–warm and hot gases. While our values for the hot |$\eta _{\rm E}$| lie within a similar range to Kim et al. (2020b) and Steinwandel et al. (2024) (⁠|$\eta _{\rm E} \approx 10^{-2}{-}10^{0}$|⁠), we find a strong negative slope in our results, unlike the positive slope from these previous works.

Bieri et al. (2023) utilize the simulating gGNs through ISM with non-equilibrium effects framework, simulating the ISM of disc galaxies at high spatial resolution, similar to the simulations previously discussed. With their implementation, Bieri et al. (2023) find mass-loading factors between 3 and 10 when measured at a distance of 5  kpc from the galaxy. These measurements are most comparable to our m12 galaxies at z|$\sim$|0, and we find that our measurements match at 750 pc scales (Fig. 2), though our galaxy-averaged measurements are slightly lower than this value (see Fig. 3). Bieri et al. (2023) also find that hot outflows are more likely to leave the galaxy than the cold and warm gas, which aligns with our velocity measurements in Fig. 6, where we show that imposing a velocity cut has no effect on the supersonic hot outflows.

The specific implementation and prescriptions of various physics can have a vast impact on outflow measurements, star cluster properties, and the effect of SNe on the ISM. For example, Gutcke et al. (2021) stress the sensitivity of the ISM to spatial energy deposition, noting that it can affect outflow behaviour and even reduce |$\eta _{\rm M}$| by a factor of 2–3 in some cases.

Recent studies also suggest that including cosmic rays can have a significant effect on measured galaxy outflows (Ruszkowski, Yang & Zweibel 2017; Hopkins et al. 2021; Chan et al. 2022). For example, Ruszkowski et al. (2017) and Hopkins et al. (2021) both find that not including cosmic rays in L|$_\star$| galaxies around can cause most outflows to simply be recycling material, while including cosmic rays leads to stellar feedback driving stronger winds. Chan et al. (2022) find a similar result, focusing specifically on the context of superbubbles in low-redshift L|$_\star$| galaxies. Rathjen et al. (2023) also suggest that cosmic rays may also dramatically affect the phase structure of outflows, as excluding them results in winds being single-phase (hot) at 2 kpc, in contrast to having a distinct multiphase structure (cold, warm, and hot) with the inclusion of cosmic rays. Furthermore, They also measure mass loading factors to be near unity with cosmic rays, similar to our z|$\sim$|0 measurements.

Runaway stars may also impact mass loading factors. Andersson, Agertz & Renaud (2020) show that, for outflows measured at 5 kpc from the galaxy in their MW–mass hydrodynamical simulations, including runaway stars may boost |$\eta _{\rm M}$| by nearly an order of magnitude. Furthermore, they find that such stars are able to transfer energy into the halo, resulting in a heating of the CGM.

Furthermore, Smith, Sijacki & Shen (2018) note that in their isolated disc galaxy setups, delayed cooling results in feedback that destroys the disc. While this feedback is matches observations in that it is effective at suppressing star formation and agrees with the Kennicutt–Schmidt relation, but they also find that only their highest resolution simulation produces large enough |$\eta _{\rm M}$|⁠. The authors specifically draw attentionto the fact that this could be because they are oversimplifying the evolution of SN remnants, need to include other stellar feedback processes, or that it could be due to their star formation prescription.

There are several reasons why we could be seeing such higher energy-loading factors in FIRE-2. One aspect could be the particular numerical implementation of SN of stellar feedback, combined with the debated burstiness of the simulated galaxies (Sparre et al. 2017; Faucher-Giguère 2018; Orr et al. 2018, 2020, 2021; Stern et al. 2021; Gurvich et al. 2023; Hopkins et al. 2023; Sun et al. 2023). The simulations by Kim et al. (2020a) and Steinwandel et al. (2024) also lack significant features including a realistic cosmological background and a pre-existing CGM.

Due to the aforementioned caveats, we stress that care should be taken when comparing our loading factors with other simulations. There clearly exists a significant need to determine the extent to which physics (e.g. cooling, star formation, turbulence, and cosmic rays), simulation parameters and methods (e.g. cosmological background, simulation scale, existing CGM, and resolution), and outflow calculations (e.g. shells versus surfaces, velocity cuts, and differing heights above the galactic centre) affect loading factors and other observed properties of outflows driven by stellar feedback.

4.1.1 Within the FIRE simulations

Muratov et al. (2015), using the original FIRE suite, found that galactic wind efficiency was dependent on halo mass, resulting in SMC/LMC-mass galaxies being extremely efficient at shutting off star formation through the expulsion of material in galactic winds. These results are aligned with what we find in the FIRE-2 galaxies, where more massive galaxies (especially those at high redshift) have less efficient galactic outflows compared to star formation on global scales (Bassini et al. 2023). FIRE low-mass galaxies maintain their bursty star formation and wind recycling up through local times (Muratov et al. 2015; Anglés-Alcázar et al. 2017; Hayward & Hopkins 2017), agreeing with the results we have found here.

Hopkins et al. (2023) provided a framework for understanding breathing modes of star formation/galactic winds: bursty star formation transitions to smoother star formation when the escape velocity becomes large enough that stellar feedback does not drive the escape of cold gas mass-loaded winds (i.e. the gravitational potential becomes deep enough to confine SN ejecta). Stellar feedback then is largely confined to the disc and star formation becomes more steady. Additionally, Stern et al. (2021) and Hafen et al. (2022) showed that as the galaxy halo grows, the inner CGM virializes, becoming hot and promoting sub-sonic gas motion, resulting in outflows that are less likely to travel significant distances from their host galaxy.

Another work using the FIRE-2 simulations, Pandya et al. (2021) found that the hot phase of gas carries both more mass and energy for MW-mass galaxies. This stands in contrast with other simulation predictions that the cold phase carries most of the mass and the hot phase carries most of the energy (Kim et al. 2020b; Fielding & Bryan 2022). In our work, we find that on the local scale, the cooler phase dominates the mass for the m12 galaxies, while the hot phase dominates the energy at high z, while the cold–warm phase carries more energy at low z (see Table 2). We see the same relationship for galaxy-averaged quantities. However, we note that we define outflows differently from Pandya et al. (2021), which use spherical shells further out at 0.1|$\rm {\it R}_{vir}$|⁠, with different velocity cuts and averaged their loadings over Gyr-time-scales

Even within the FIRE-2 simulation suite, we note that resolution can affect galaxy properties analysed here, including gas fractions, stellar masses, and star formation, as physical processes that affect such quantities can be resolution dependent

Hopkins et al. (2018) conduct an in-depth analysis into the effects of resolution in Section 4, specifically presenting tests of mass resolution in their figs 8–14. In particular, they note that in MW-mass galaxies, including one of the m12 galaxies used here, increasing resolution can shift star formation to later times as the hot gas and star formation in progenitor galaxies are resolved. The SFR is therefore affected due to the changes in the gas reservoirs of the CGM. Similarly, Hopkins et al. (2020) find that lower resolution in the m12 galaxies can produce larger stellar masses, enhanced circular velocities, and higher central baryonic densities.

Garrison-Kimmel et al. (2019) also find that mass resolution can significantly affect the star formation of low-mass galaxies, such as our m11s. For example, when run with a higher resolution of |$m_{\rm b, min}=880 \rm\, {\rm M}_{\odot}$|⁠, m11h at z|$\sim$|0 ends up with |$M_{\star } \approx 10^{8.1} \rm\, {\rm M}_{\odot}$| and an irregular morphology, whereas at |$m_{\rm b, min}=7100 \rm\, {\rm M}_{\odot}$| (as used in this paper), the final stellar mass is |$M_{\star }\approx 10^{9.6} \rm\, {\rm M}_{\odot}$| with a rotationally supported disc (El-Badry et al. 2018; Porter et al. 2022). Because the higher resolution run of m11h has a lower stellar mass at z|$\sim$|0, it also has a SFR that is nearly 2.5 dex lower than that of the |$m_{\rm b, min}=7100 \rm\, {\rm M}_{\odot}$| run.

Because changes in mass resolution for some of the FIRE-2 zoom-in simulations studied here can result in different final galaxies and/or global galaxy properties, such as in the case of m11h, or m12i (see Hopkins et al. 2018; Gurvich et al. 2020), we rationalize our decision to only compare results across the |$m_{\rm b, min}=7100 \rm\, {\rm M}_{\odot}$| m11 and m12 simulations.

4.2 Comparison to recent observations

In Fig. 3, we compare directly with deep H|$\alpha$| observations of outflows in six starbursting galaxies by McQuinn et al. (2019). They suggest that true outflowing feedback-driven winds are preferentially identified in low-mass galaxies with centrally concentrated star formation, similar to the findings by Fielding et al. (2018) and Orr et al. (2022a), who predict feedback from high gas surface density central regions produce winds more likely to escape the galaxy. Indeed, we identified such winds in the SMC-mass galaxy m11d, that does host centrally clustered star formation (see the centre panels of Fig. 7).

We also compare the |$\sim$|kpc-scale momentum flux in the outflows of our MW-mass galaxies with similarly spatially resolved observations of IRAS08, a starbursting disc galaxy, by (Reichardt Chu et al. 2022), as described in Section 3.2. Though we measure outflowing gas in all of our m12s, it is clear to see that the FIRE-2 disc galaxies generally do not form massive nuclear star clusters capable of hosting fast SN-driven winds, the likes of which are found in IRAS08 or, for example, M82 (Leroy et al. 2015). Most of the spatially resolved ‘winds’ in the MW-mass galaxies (⁠|$M_\star \sim 10^{10}$| M|$\odot$|⁠) originate from the galactic outskirts (see Fig. 9), though where there is overlap with the distribution of momentum fluxes seen in IRAS08, we find general agreement in their magnitude.

Feedback models for regulating the ISM have highlighted the trade-off between driving turbulence in the ISM and galactic outflows (Ostriker et al. 2010; Faucher-Giguère et al. 2013; Hayward & Hopkins 2017; Krumholz et al. 2018; Orr et al. 2022b). Our results, showing a lack of strong wind driving from the central regions, are in agreement with the superbubble breakout prediction of Orr et al. (2022b).

Recent studies from ALMA and JWST have also directly helped constrain theoretical predictions of superbubble feedback. Work with JWST in NGC628 by Mayya et al. (2023) uncovered the temporally extended star formation episode driving a superbubble shell (the |$\sim$|kpc hole in the outer disc of that galaxy), in-line with our results in Fig. 8 finding that stellar feedback persists (i.e. some evidence of ‘positive feedback’) during the galaxy’s identified superbubble episode. This result is likely capturing the continued stellar feedback that contributes to driving the bubble itself, capturing an ongoing pattern of stellar feedback, though the observations of NGC 628 may indeed be capturing triggered star formation.

One important aspect we have not taken into account in this paper is the importance of viewing angles when measuring galactic outflows. Bordoloi et al. (2014) specifically measure the cool outflow rates and equivalent widths from Mg ii absorption lines. They find that the equivalent width of these absorption lines depends strongly on the star formation rate surface density, and that the equivalent width of this line in disc galaxies is largest for face-on galaxies, suggesting that outflows have a bipolar geometry and velocities are greater in face-on galaxies than edge-on. They also find that irregular morphologies, such as the low-mass galaxies in our sample, have outflow equivalent widths and velocities similar to the face-on discs. The importance of viewing angle is also supported by recent JWST NIRCam and NIRSpec observations at higher redshifts (z|$\sim$| 6–9), where Zhang et al. (2024) suggest that one possibility for galaxies with (or without) extended gas emission or broad lines can be explained by different viewing angles with outflows. While investigating this topic is outside the scope of this paper, simulations testing different viewing angles and resulting outflow measurements/line widths could prove a valuable future comparison.

5 SUMMARY AND CONCLUSIONS

In this paper, we presented an analysis of supernova feedback-driven outflows in 11 FIRE-2 galaxies, six of which are MW-mass progenitors and five of which are SMC/LMC-mass progenitors across a redshift range from |$z \approx 3{-}0$|⁠. At higher redshifts (⁠|$z\ge 1$|⁠), we calculate outflows through a surface at a distance |$\rm \pm 0.05{\it R}_{vir}$| above/below the galaxy, while at |$z=0$| we calculate outflows at two galaxy-averaged gas scale heights above/below the galaxy. We note that both choices have similar values (see Appendix  A).

Comparing our findings from the FIRE suite with other simulations and recent observations, we found good agreement with observed and simulated (Kim et al. 2020b; Bieri et al. 2023; Steinwandel et al. 2024) mass loadings of winds, but elevated energy loadings compared to the TIGRESS and Steinwandel et al. (2024) simulations. However, we note that these simulations differ in fundamental ways that may affect measurements, including calculation of star formation rates and cosmological context (isolated galaxies and ISM patches lacking a CGM), as detailed in Section 4.1. We picked out several remarkable instances of galaxy-wide superbubbles occurring in simulations, and characterized the behaviour of their resulting outflows.

Our main findings can be summarized by the following points:

• On pixel (spatial) scales of 750 pc, we find negatively sloped power laws in both spatially resolved mass loadings |$\eta _{\rm M}$| (Fig. 2) and energy loadings |$\eta _{\rm E}$| (Fig. 4) versus. star formation rate surface density |$\dot{\Sigma }_\star$|⁠, with the steepest slopes found at |$z\approx 0$|⁠. Our values for the cold–warm gas |$\eta _{\rm M}$| are comparable with observations from McQuinn et al. (2019) and simulations by Steinwandel et al. (2024) and (Kim et al. 2020b). We predict higher hot gas mass loadings across redshift than either the TIGRESS simulations or those of Steinwandel et al. (2024).

• On 750 pc scales (Figs 2 and 4), the MW-mass galaxy progenitors (m12s) and SMC-mass galaxies (m11s) have nearly identical mass-loadings |$\eta _{\rm M}$| as a function of star formation rate surface density. The MW-mass galaxies do have higher energy loadings |$\eta _{\rm E}$| on these scales, which may be due to the halo potentials; the galaxy must launch much more energetic outflows to reach the same defined height (2H, 0.05|$\rm {\it R}_{vir}$|⁠) as the smaller m11s.

• On galaxy-averaged scales (Figs 3 and 5), the SMC-mass m11 galaxies have higher mass-loading values than the MW-mass progenitors (m12s), especially at z|$\sim$|0. The mass loading is also dominated by the cold–warm phase at all points in time. The m11s also appear to have higher energy loading at z|$\sim$|0, which may simply be because the m11 galaxies, though launching similar mass/energy-loaded winds locally, have lower global SFRs. At high redshifts for both galaxy groups, the energy loading is dominated by the hot gas, while at low z the energy distribution is closer to equipartition.

• These simulated galaxies all exhibit bursty star formation modes prior to the formation of discs. Following starburst events, in nearly all cases, we see significant outflowing mass and energy flux originating from clearly identifiable superbubble events.

• Most high-redshift snapshots of all galaxies in our analysis have galaxy-averaged ISM properties that correspond to the ‘Powered Breakout’ regime from Orr et al. (2022a), with local values that are comparable to observations from DUVET (Reichardt Chu et al. 2022), suggestive of the idea that most if not all of these simulations host star-forming regions capable of driving outflows.

• Variation in measured outflow rates, in both |$\eta _{\rm M}$| and |$\eta _{\rm E}$|⁠, can vary significantly (an order of magnitude of more) depending on where outflows are measured (see Fig. A1), and thermal motions from the inner CGM may play a significant role in these measurements.

Quantifying the effects of clustered stellar feedback on galaxies of all types across cosmic time is essential to grasping what modes of star formation are capable of driving the strong galactic outflows seen in the universe. We know these outflows have direct consequences on galaxy formation: affecting accretion, star formation, the mass–metallicity relation, etc., and a firm understanding of these processes directly constrains our theories of galaxy formation and evolution. Future theoretical work understanding galaxy-scale superbubbles and feedback-driven winds will become only more necessary as observational constraints on mass and energy loadings at all redshifts tighten, and more detailed kinematics of superbubbles are measured.

ACKNOWLEDGEMENTS

We would like to thank Ulrich Steinwandel, Drummond Fielding, and the anonymous referee for useful conversations and comments that greatly improved the manuscript. We ran simulations using: XSEDE, supported by National Science Foundation (NSF) grant ACI-1548562; Blue Waters, supported by the NSF; Frontera allocations AST21010 and AST20016, supported by the NSF and Texas Advanced Computing Center (TACC); Pleiades, via the National Aeronautics and Space Administration (NASA) High-End Computing (HEC) program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. LEP is thankful to the Flatiron Institute for continued financial support and resources. The Flatiron Institute is supported by the Simons Foundation. BB is grateful for generous support by the David and Lucile Packard Foundation and Alfred P. Sloan Foundation. AW received support from NSF via CAREER award AST-2045928 and grant AST-2107772; NASA ATP grant 80NSSC20K0513. C-AF-G was supported by NSF through grants AST-2108230, AST-2307327, and CAREER award AST-1652522; by NASA through grants 17-ATP17-0067 and 21-ATP21-0036; and by the Space Telescope Science Institute (STScI) through grants HST-GO-16730.016-A and JWST-AR-03252.001-A.

DATA AVAILABILITY

The data supporting the plots within this article are available on reasonable request to the corresponding author. A public version of the gizmo code is available at http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html. Additional data including simulation snapshots, initial conditions, and derived data products are available at https://fire.northwestern.edu/data/. The FIRE-2 simulations are publicly available (Wetzel et al. 2023) at http://flathub.flatironinstitute.org/fire.

Footnotes

REFERENCES

APPENDIX A: DEPENDENCE OF OUTFLOW FLUXES ON HEIGHT/DISTANCE OF MEASURED FLUX SURFACE

In Fig. A1, we show how the various loading factors change if we choose to place the surface through which outflows are defined at different heights above the galaxy at z|$\sim$|1–3. We do not distinguish between galaxy type or impose a cut on velocity at the measured surface, instead focusing solely on the differences in height at which outflows are measured.

We can see that there is little difference in flux between any of the heights (⁠|$\rm H$|⁠, |$\rm 2H$|⁠, |$\rm 0.05{\it R}_{vir}$|⁠, and |$\rm 0.1{\it R}_{vir}$|⁠) for the hot gas mass-loading factor (upper right panel). Unsurprisingly, Fig. A1 shows a difference of up to an order of magnitude between the lowest surface height, |$\rm H$|⁠, and the highest surface height, |$\rm 0.1{\it R}_{vir}$|⁠, in a consistent manner: There is more flux through surfaces closer to the galaxy mid-plane. We include this comparison here to show to what extent our fiducial chosen heights affects our results.

While we have shown that the relationship between loading factors and the star formation rate surface density are affected by the chosen surface height, we stress that we frequently analyse how the gas and energy flux differ over a period of time where incidents of strong stellar feedback have been identified (see e.g. Fig. 8). It becomes important in these cases that the surface where we define outflows is not only at an appropriate height compared to the galaxy as a whole, but that this surface height is consistent. Otherwise, our measurements of outflow properties may not reflect the true effect of feedback. For example, if we define the measured outflow surface at the gas scale height, when a single massive superbubble causes gas to be disrupted on a galaxy-wide scale, the gas scale height can dramatically change from snapshot to snapshot. If we measure outflow properties through the surface one or two scale heights above the galaxy during one of these events, the outflow properties might just reflect the growing/expanding ‘cold cap’ of the superbubble as the galaxy quite literally breaks apart. In contrast, measurements at fractions of the virial radius are much more stable against large-scale disruptions as the virial radius grows smoothly with cosmic time.

APPENDIX B: ENERGY FLUXES AND STAR FORMATION RATES FOR INDIVIDUALLY IDENTIFIED SUPERBUBBLE EVENTS

Fig. B1 shows the evolution in the energy flux of the hot gas component as presented in Fig. 8, but for all identified superbubbles in all individual galaxies of this study.

These events occur in all four redshift bins (⁠|$z\sim 3{-}0$|⁠) and both galaxy types of this study. In addition, we note that the same relationship is observed in every single identified superbubble: The star formation rate reaches a peak (typically for the entire time plotted) during the superbubble snapshots, then quickly drops. This drop in star formation is between 0.5 and 1 dex, with the exception of m12b, which experiences a smaller drop of 0.25–0.5 dex.

Table B1 shows the relative percentages of outflows from each identified superbubble, where the mass and energy flux from each superbubble instance (one snapshot each for before, during, and after) is divided by the total mass and energy flux from that galaxy over the specific redshift range. If superbubbles did not contribute significantly to the total flux over this time period, we might expect each snapshot to contribute about 10 per cent of the total. Therefore, we can consider the contributions from superbubbles significant compared to the total if the superbubble snapshots contribute more than 10 per cent of the flux times the number of snapshots identified to have a possible superbubble (approximately three per superbubble episode). In this case, this would apply to m11d at z|$\sim$|3, m11d at z|$\sim$|2, m11d at z|$\sim$|0, m11e at z|$\sim$|2, m12r at z|$\sim$|3, and m12r at z|$\sim$|1.