Code Comparison in Galaxy-scale Simulations with Resolved Supernova Feedback: Lagrangian versus Eulerian Methods

We use four different codes for gravity and hydrodynamics: one Eulerian code (Ramses), one moving mesh code (Arepo), and two Lagrangian codes (Gizmo and Gadget-3), which we briefly describe as follows.

Gadget-3 (hereafter Gadget for short) (Springel 2005) is a particle-based code. Gravity is solved using a "tree code" (Barnes & Hut 1986) where short-range forces are calculated pairwise while long-range forces are approximated by the center of mass of a tree node. Since it is designed for collisionless N-body problems, gravitational forces at small distances are reduced by softening, and different softening lengths can be used for different components, such as gas, stars, and dark matter. Hydrodynamics is solved by the smoothed particle hydrodynamics (SPH) method (Lucy 1977). We adopt the SPHGal implementation in Hu et al. (2014) that includes several improvements over traditional SPH methods; these improvements include the pressure-energy formulation (Hopkins 2015), the Wendlend kernel (Dehnen & Aly 2012), variable artificial viscosity and conduction (Price 2008; Cullen & Dehnen 2010), and the time-step limiter (Durier & Dalla Vecchia 2012). Both gravity and hydrodynamics solvers are Lagrangian in the sense that particles follow the motions of the gas, stars, and dark matter. The initial mass carried by each particle is conserved as there are no mass fluxes between particles.

Gizmo (Hopkins 2015) is a multimethod code built on the code Gadget, featuring the meshless Godunov method (Gaburov & Nitadori 2011) for hydrodynamics. We adopt the meshless finite-mass (MFM) method proposed by Hopkins (2015) in this work. In MFM, the simulation domain is divided into overlapping cells defined by the particle distribution and a kernel function with a finite support radius. The Riemann problem is solved at the interfaces between cells. The cells move with the gas so that the mass fluxes between cells are assumed to be zero. Gizmo adopts the same gravity solver as Gadget and is a Lagrangian code when the MFM solver is used.

Arepo (Springel 2010; Pakmor et al. 2016; Weinberger et al. 2020) uses a second-order accurate finite-volume scheme on an unstructured, moving mesh. The simulation domain is discretized by a Voronoi tessellation, and a Riemann problem is solved at the interfaces between cells to compute fluxes. The discrete set of mesh-generating points that define the tessellation is allowed to move with a velocity close to that of the local fluid (with small corrections to maintain cell regularity). The method is therefore pseudo-Lagrangian, as the moving mesh tends to minimize mass fluxes between cells, with the result that they roughly maintain a constant mass. However, these mass fluxes are nonzero (in contrast to Gadget and Gizmo). In the standard usage of the code, adopted here, a refinement (de-refinement) scheme is used to split (merge) cells in order to enforce a constant mass resolution within a factor of 2. Arepo adopts a gravity solver similar to those employed in Gadget and Gizmo.

Ramses (Teyssier 2002) is an Eulerian code with octree-based adaptive mesh refinement (AMR). The N-body solver is based on the adaptive-particle-mesh method. The Poisson equation is solved using the multigrid method with Dirichlet boundary conditions at level boundaries (Guillet & Teyssier 2011). Hydrodynamics is solved using the MUSCL scheme, a second-order finite-volume Godunov method, and the HLLC Riemann solver. For Ramses, we used here a quasi-Lagrangian refinement criterion based on a target mass common to all codes in this work. We use a box size of 450 kpc with a minimum level of refinement ${{\ell }}_{\min }=8$ and a maximum level of refinement ${{\ell }}_{\max }=16$ (${{\ell }}_{\max }=17$) for the low (high) resolution run, resulting in a maximum spatial resolution of 7 pc (resp. 3.5 pc).

The mesh in Eulerian codes is static by definition. The quasi-Lagrangian AMR leads to an adaptive spatial resolution down to a minimum cell size. In contrast, in both moving mesh and Lagrangian codes, the resolution elements follow the motions of the gas, leading to a smoothly (and, in principle, infinitely) adaptive spatial resolution with a constant or near-constant mass resolution. Despite being pseudo-Lagrangian, moving mesh codes share many properties with Lagrangian codes. For convenience, we will refer to Gadget, Gizmo, and Arepo as "Lagrangian codes" hereafter.

The initial conditions are generated with the MakeDiskGalaxy code developed in Springel et al. (2005), which consists of a rotating disk galaxy embedded in a dark matter halo. The halo has a virial radius R_vir = 45 kpc and a virial mass M_vir = 10¹⁰ M_⊙, and it follows a Hernquist profile, which matches an NFW (Navarro et al. 1997) profile at small radii with the concentration parameter c = 15 and the spin parameter λ = 0.035. The baryonic mass fraction is 0.8%, out of which 10⁷ M_⊙ is in the stellar disk and 7 × 10⁷ M_⊙ in the gaseous disk. Both the stellar and gaseous disks follow an exponential profile a with scale length of 1 kpc, which makes the central gas surface density Σ_gas ∼ 10 M_⊙ pc⁻². The stellar disk has a scale height of 1 kpc, while the vertical profile of the gaseous disk is such that it is in vertical hydrostatic equilibrium. The initial gas temperature is set to 10⁴ K. The hydrogen mass fraction is X_H = 0.76. Our initial conditions resemble the Wolf–Lundmark–Melotte (WLM) galaxy, a nearby star-forming dwarf irregular galaxy that has multiwavelength observations (Leaman et al. 2012; Rubio et al. 2015; Mondal et al. 2018). The simulation is run for 1 Gyr for the Lagrangian codes and 0.5 Gyr for Ramses due to computational cost.

For Ramses and Arepo, it is necessary to set up a minimum density for numerical purposes, which is n = 10⁻⁷ cm⁻³ uniformly distributed in the background. In contrast, there is no background gas in the halo for Gizmo and Gadget.

For comparison purposes, we should, ideally, adopt the same resolution in both Lagrangian and Eulerian codes. However, Lagrangian codes have a fixed mass resolution and an adaptive spatial resolution while Eulerian codes are the opposite. Due to this intrinsic difference, there is no unique way of choosing the same resolution for both methods. In this work, we adopt a spatial resolution of Δx = 7 pc in the low-resolution model for the Eulerian code Ramses. At this resolution, the thermal Jeans length can be resolved up to n ∼ 100 cm⁻³ assuming T = 30 K. Therefore, for the Lagrangian codes, we choose the gas particle mass m_g such that the smoothing length is comparable to Δx at n = 100 cm⁻³, which translates to m_g = 100 M_⊙ in the low-resolution models. The gravitational softening for the baryons (i.e., gas, preexisting stars, and the stars formed in the simulation) is set to be h_g = 10 pc. The dark matter halo is resolved with a much coarser resolution with the particle mass m_dm = 10⁴ M_⊙ and the gravitational softening h_dm = 200 pc. In the high-resolution models, we simply decrease all the spatial resolutions by a factor of 2 and the mass resolutions by 8. Table 1 summarizes the resolution-related parameters.

We use the public Grackle library (Smith et al. 2017) ¹² for radiative cooling, adopting its equilibrium cooling table without solving a chemistry network. The gas metallicity is Z = 0.002 (0.1 Z_⊙) and is constant throughout the simulation. In addition, we include heating from the photoelectric effect. Following Bakes & Tielens (1994) and Wolfire et al. (2003), the photoelectric heating rate can be expressed as

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{PE}}=1.3\times {10}^{-24}n{\epsilon }_{\mathrm{PE}}{G}_{0}{Z}_{{\rm{d}}}^{{\prime} }\,\mathrm{erg}\,{\mathrm{cm}}^{-3}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 1 }$$

where n is the hydrogen number density, _PE is the photoelectric heating efficiency, G₀ is the radiation strength in units of the Habing field (Habing 1968), and ${Z}_{{\rm{d}}}^{{\prime} }$ is the dust-to-gas mass ratio relative to the Milky Way value (∼1%). We adopt a constant _PE of 0.05 (see Appendix A) and ${Z}_{{\rm{d}}}^{{\prime} }=0.1$, assuming a linear Z–${Z}_{{\rm{d}}}^{{\prime} }$ relationship. ¹³ We assume G₀ scales linearly with the SFR surface density (Σ_SFR), normalized to the solar-neighborhood values of G₀ = 1.7 and Σ_SFR = 2.4 × 10⁻³ M_⊙ yr⁻¹ kpc⁻². Since the SFR is unknown prior to the simulation, we use low-resolution simulations to empirically determine the photoelectric heating rate as a function of the galactocentric radius R:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{PE}}={{\rm{\Gamma }}}_{\mathrm{PE},0}\exp \Space{0ex}{0.25ex}{0ex}(\displaystyle \frac{{R}_{{\rm{c}}}-\max (R,{R}_{{\rm{c}}})}{{R}_{{\rm{c}}}}\Space{0ex}{0.25ex}{0ex}),\end{eqnarray} \tag{ 2 }$$

where Γ_PE,0 = 2.6 × 10⁻²⁷ erg s⁻¹ and R_c = 0.4 kpc. Our low-resolution experiments suggest that the results are insensitive to modest variations of Γ_PE (a factor of 3), which is consistent with Hu et al. (2017) and Smith et al. (2021).

We adopt the stochastic star formation recipe commonly used in simulations of galaxy formation. The local SFR is ${\dot{\rho }}_{* }={\epsilon }_{\mathrm{SF}}{\rho }_{\mathrm{gas}}/{t}_{\mathrm{ff}}$, where ρ_gas is the gas density, _SF is the star formation efficiency, and ${t}_{\mathrm{ff}}\equiv \sqrt{3\pi /(32G{\rho }_{\mathrm{gas}})}$ is the free-fall time, and G is the gravitational constant. We adopt _SF = 1% as our fiducial choice, but we also explore _SF = 100% in some models. Instead of adopting a density threshold for star formation, we allow star formation only to occur when the thermal Jeans length ${L}_{{\rm{J}}}=\sqrt{\pi {c}_{s}^{2}/(G{\rho }_{\mathrm{gas}})}\lt {L}_{{\rm{J}},0}$, where c_s is the sound speed and L_J,0 is a predefined star formation threshold. We adopt L_J,0 = Δx such that gas is eligible for star formation only when the Jeans length becomes unresolved as we can no longer follow the gravitational collapse faithfully (Truelove et al. 1997). For Lagrangian codes, this star formation criterion roughly corresponds to the density where the Jeans mass becomes unresolved (Bate & Burkert 1997).

Star particles are stochastically created based on the local SFR. In the Lagrangian codes, the gas resolution element is converted into a star particle that inherits the mass of its parent. This means that all star particles in Gadget and Gizmo have a mass of m_g exactly, while star particles in Arepo are guaranteed to have a mass within a factor of 2 of m_g. In Ramses, a star particle with a mass of m_g is spawned by removing a gas mass of m_g from the parent cell.

We include feedback from core-collapse SNe. For a typical stellar population, there is about one SN progenitor in every 100 M_⊙ stellar mass. As our star particle has a mass of m_* ≤ 100 M_⊙, there is at most one SN progenitor in each star particle statistically. Therefore, we stochastically sample massive stars such that each star particle has a probability of m_*/(100 M_⊙) of being selected as an SN progenitor with an SN delay time of ${t}_{\mathrm{SN}}=10\,\mathrm{Myr}$. Each SN injects 10⁵¹ erg of thermal energy into the surrounding ISM. For Ramses and Arepo, the energy is injected into the cell where the star is located. For Gizmo and Gadget, multiple cells can overlap with the same star simultaneously. For simplicity, we inject energy into the nearest eight particles in a kernel-weighted fashion. For each SN event, we record not only the ambient gas density and temperature but also the location and time such that we can perform a cluster analysis in Section 3.4.

The models we run are summarized in Table 2. For the low-resolution case, we only run our fiducial model where _SF = 1% and ${t}_{\mathrm{SN}}=10\,\mathrm{Myr}$. For the high-resolution case, we run our fiducial model, a model with ${t}_{\mathrm{SN}}=0$, and a model with _SF = 100%. For Ramses, we run an additional model with ${t}_{\mathrm{SN}}=50\,\mathrm{Myr}$. Finally, we run a model for the Lagrangian codes (Gizmo, Arepo, and Gadget), where we adopt ${t}_{\mathrm{SN}}=0$ up to t = 250 Myr and ${t}_{\mathrm{SN}}=10\,\mathrm{Myr}$ afterward. As we will demonstrate, this particular setup turns out to be our favored fiducial model for Lagrangian codes as it prevents the artificial blowout of the disk during the initial phase. We do not have the Gadget runs for all the models (i.e., there is no gadget and gadget_sfe100).

Table 2. Overview of Simulation Models

Model Name	Code	Resolution	SF	${t}_{\mathrm{SN}}$ (Myr)	Additional Comments
gizmo_lr	Gizmo	100 M⊙	1%	10
arepo_lr	Arepo	100 M⊙	1%	10
gadget_lr	Gadget	100 M⊙	1%	10
ramses_lr	Ramses	7 pc	1%	10
gizmo	Gizmo	12.5 M⊙	1%	10
arepo	Arepo	12.5 M⊙	1%	10
ramses	Ramses	3.5 pc	1%	10
gizmo_tsn0	Gizmo	12.5 M⊙	1%	0
arepo_tsn0	Arepo	12.5 M⊙	1%	0
gadget_tsn0	Gadget	12.5 M⊙	1%	0
ramses_tsn0	Ramses	3.5 pc	1%	0
ramses_tsn50	Ramses	3.5 pc	1%	50
gizmo*	Gizmo	12.5 M⊙	1%	10	${t}_{\mathrm{SN}}=0$ until t = 0.25 Gyr
arepo*	Arepo	12.5 M⊙	1%	10	${t}_{\mathrm{SN}}=0$ until t = 0.25 Gyr
gadget*	Gadget	12.5 M⊙	1%	10	${t}_{\mathrm{SN}}=0$ until t = 0.25 Gyr
gizmo_sfe100	Gizmo	12.5 M⊙	100%	10
arepo_sfe100	Arepo	12.5 M⊙	100%	10
ramses_sfe100	Ramses	3.5 pc	100%	10

Note. Models with an asterisk are run with ${t}_{\mathrm{SN}}=0$ for the first 0.25 Gyr to prevent the initial artificial blowout of the gaseous disk. _SF: star formation efficiency. ${t}_{\mathrm{SN}}$: SN delay time.

Download table as:  ASCIITypeset image

Figure 1 shows the star formation rate (SFR) as a function of time for the low- and high-resolution models in the top and bottom panels, respectively. The average SFRs in the low-resolution models are in broad agreement in all codes. However, the ramses_lr model is notably less bursty (i.e., the SFR varies less with time) than the Lagrangian models. In addition, there is an initial starburst in the Lagrangian models at t ∼ 100 Myr that is absent in the ramses_lr model. In the high-resolution models (bottom panel), the SFRs in the Lagrangian models become even burstier. Indeed, the initial starburst leads to extremely energetic SN feedback, which blows out the entire gaseous disk. Although a quasi-steady state is eventually established after some of the blown-out gas falls back, the overall gas surface density becomes significantly lower than that in the initial conditions, leading to a reduced time-averaged SFR. The blowout is, however, a somewhat artificial consequence of our initial conditions and numerical setup. Before the first star formation occurs, gas cools and collapses globally into an artificially thin disk. Once the gas becomes dense enough to form stars, the resulting coherent SNe can easily break out of the thin disk, resulting in a catastrophic blowout. This is a well-known phenomenon in simulations of isolated galaxies as well as in ISM patches, and several methods have been adopted to mitigate the initial blowout, such as turbulent driving (e.g., Walch et al. 2015; Kim & Ostriker 2017) or random SN driving (e.g., Hu et al. 2017; Smith et al. 2021). On the other hand, the ramses model still shows a nonbursty SFR similar to its low-resolution counterpart and there is no initial blowout of the disk, in sharp contrast to the Lagrangian models.

Zoom In Zoom Out Reset image size

An important clue in understanding this difference comes from the observation that the SN delay time has a significant effect on the burstiness of star formation in the Lagrangian codes. The top panel of Figure 2 shows the SFR time evolution for models with instantaneous SN feedback (${t}_{\mathrm{SN}}$ = 0). With this change, the SFRs in gizmo_tsn0, ramses_tsn0, and gadget_tsn0 all show excellent agreement with each other. The SFR in arepo_tsn0 is lower by a factor of 2, and it is unclear what causes the slight difference. More importantly, the SFRs in all codes are strikingly constant with little temporal fluctuation. In addition, the initial burst is absent and a quasi-steady state is rapidly established without the blowout of the disk.

Zoom In Zoom Out Reset image size

An interesting implication is that instantaneous SN feedback can be a simple alternative method to mitigate the artificial initial starburst. Motivated by this, we rerun our Lagrangian codes with ${t}_{\mathrm{SN}}$ = 0 for t < 250 Myr and ${t}_{\mathrm{SN}}$ = 10 Myr for t ≥ 250 Myr. We do not do so with Ramses as it shows no initial starburst. The results are shown in the middle panel of Figure 2. The Lagrangian models show reasonable agreement with each other. Although there is still a burst of SFR at t ≳ 250 Myr right after we switch ${t}_{\mathrm{SN}}$ from 0 to 10 Myr, it is much weaker and does not blow out the disk. Consequently, the time-averaged SFR is not reduced and is comparable to ramses. We therefore refer to these models as our "fiducial" Lagrangian models. However, even if the artifacts of the initial conditions have been greatly reduced, the SFRs for the Lagrangian codes are still burstier than ramses. This suggests that the difference in burstiness between the two methods is an intrinsic property rather than an artifact of the initial conditions.

The observed SFR in the WLM galaxy from Hunter et al. (2010) is 1.7 × 10⁻³ M_⊙ yr⁻¹ from Hα and 6.3 × 10⁻³ M_⊙ yr⁻¹ from far-ultraviolet (FUV). The SFRs in our fiducial models are broadly in good agreement with observations, which is reassuring.

The bottom panel of Figure 2 shows models with _SF= 100% as well as a Ramses model with ${t}_{\mathrm{SN}}$ = 50 Myr. Both gizmo_sfe100 and arepo_sfe100 show very bursty SFRs similar to the corresponding models gizmo and arepo where _SF = 1%. On the other hand, both ramses_sfe100 and ramses_tsn50 are bursty, in contrast to ramses. In other words, the burstiness in Eulerian codes can be enhanced by increasing either _SF or ${t}_{\mathrm{SN}}$. We will explore the reason for this in later sections.

The time-averaged SFRs of different models are summarized in Table 3.

Figure 3 shows the face-on maps of the gas surface density in different models at t = 500 Myr. The most striking feature is that large SN-driven bubbles are clearly seen in some models but are completely absent in others. In fact, the presence of SN bubbles is closely correlated with the burstiness of the SFR shown in Figures 1 and 2: Models with burstier star formation show larger SN bubbles. In particular, those that experience an initial blowout have the largest SN bubbles, and their gas surface densities are notably reduced even at t = 500 Myr, as some of the blown-out gas never falls back to the disk. On the other hand, SN bubbles are completely absent in models where ${t}_{\mathrm{SN}}$ = 0. Similarly, the standard ramses runs show very few SN bubbles, consistent with its nonbursty star formation history. SN bubbles in Ramses models can only be generated by increasing either _SF (ramses_sfe100) or ${t}_{\mathrm{SN}}$ (ramses_tsn50).

Zoom In Zoom Out Reset image size

Figure 4 shows the two-dimensional mass-weighted normalized distribution of the hydrogen number density versus temperature (the so-called "phase diagram") at t = 500 Myr. The red solid line indicates the star formation threshold where L_J = 3.5 pc. Broadly speaking, all models show similar distributions in the phase diagram despite their differences in star formation burstiness and gas morphology. This is because the thermal balance in the ISM is mainly controlled by radiative cooling. The majority of gas is in the diffuse warm gas, where n ∼ 0.1 cm⁻³ and T ∼ 10⁴ K, which is consistent with Hu et al. (2016, 2017). However, models with ${t}_{\mathrm{SN}}=0$ produce less hot (T > 10⁵ K) and diffuse gas compared to their fiducial counterparts. This is expected as the hot and diffuse gas typically exists in the SN bubbles that are largely absent in these models. The Ramses runs have more gas in the density range of n ∼ 10–100 cm⁻³ as indicated by the slightly brighter color but have very little gas above the star formation threshold.

Zoom In Zoom Out Reset image size

The environment where SNe occur provides crucial information on the efficiency of SN feedback. The left panel of Figure 5 shows the cumulative distribution of the ambient gas density where SNe occur (${n}_{\mathrm{SN}}$) for the Lagrangian models. The agreement between different Lagrangian codes is remarkable. Most SNe occur in diffuse gas where ${n}_{\mathrm{SN}}\sim {10}^{-2}\,{\mathrm{cm}}^{-3}$ in the fiducial models (gizmo*, arepo*, and gadget*) and therefore are well resolved. In contrast, most SNe occur in dense gas where ${n}_{\mathrm{SN}}\sim {10}^{2}\,{\mathrm{cm}}^{-3}$ in the instantaneous SN models (gizmo_tsn0, arepo_tsn0, and gadget_tsn0). This is expected, as SN feedback kicks in right after star formation, which by construction only occurs in cold and dense gas, ¹⁴ leading to more efficient energy loss from radiative cooling, which explains the nonburstiness and the lack of SN bubbles. Although these SNe are formally unresolved, the time-averaged SFRs in these models are still comparable to those in the fiducial models, suggesting that these "unresolved" SNe are still able to inject enough momentum to regulate star formation.

Zoom In Zoom Out Reset image size

In contrast, ramses shows a very different ${n}_{\mathrm{SN}}$ distribution compared to the Lagrangian codes with a significant fraction of SNe occurring at n > 10 cm⁻³, as shown in the right panel of Figure 5. This is part of the explanation for why ramses behaves more like the instantaneous SN models. However, there are still ∼25% of SNe with ${n}_{\mathrm{SN}}\lt 0.01\,{\mathrm{cm}}^{-3}$ in ramses, which is not significantly lower compared to gizmo*. This implies that SNe occurring at low densities is not a sufficient condition for efficient feedback, which we will discuss in more detail in Section 3.4. On the other hand, both ramses_sfe100 and ramses_tsn50 have a large fraction of low-${n}_{\mathrm{SN}}$ SNe and therefore show busty SFRs and large SN bubbles.

Having established that the burstiness and gas morphology are both related to ${n}_{\mathrm{SN}}$, we now need to understand what leads to the difference in ${n}_{\mathrm{SN}}$ in different models.

The left panel of Figure 6 shows the cumulative distribution of the gas density where star formation occurs (n_SF) for the Lagrangian codes. The agreement between different Lagrangian codes is again remarkable. With instantaneous SN feedback (dashed lines), most star formation occurs in a narrow range of densities where n ∼ 500 cm⁻³. Once the gas reaches high enough densities to form stars, SN feedback that occurs instantaneously is able to stop the gas from further collapsing. Therefore, stars always form around the star formation threshold density. In contrast, in the fiducial models with ${t}_{\mathrm{SN}}=10\,\mathrm{Myr}$ (solid lines), star formation occurs over a large range of densities, spanning more than two orders of magnitude from 500 to 10⁵ cm⁻³. In this case, gas keeps collapsing toward higher densities even after it exceeds the threshold density, as there is no countering mechanism against gravity before the first SNe occur in 10 Myr. In contrast, as shown in the right panel of Figure 6, ramses behaves very differently from its Lagrangian counterparts, showing a narrow range of n_SF very similar to the instantaneous SN models.

Zoom In Zoom Out Reset image size

We therefore conclude that Lagrangian and Eulerian codes behave very differently at densities above the star formation threshold where the Jeans length becomes unresolved: While gas in Lagrangian codes continues to collapse to much higher (unresolved) densities, gas in Eulerian codes lingers around the star formation threshold. This is because the spatial resolution in Lagrangian codes is adaptive and therefore it is much easier for the unresolved gas to collapse. In comparison, the collapse of unresolved gas in Eulerian codes is limited by the fixed spatial resolution and is thus relatively more difficult. Note that in both cases, the collapse is unresolved (and additional physical processes such as pre-SN feedback may become important) and so neither behavior is necessarily correct.

In this section, we study the clustering properties of SNe and show that SN clustering has a fundamental impact on ${n}_{\mathrm{SN}}$, SFR burstiness, and gas morphology. The effect of SN clustering is two-fold. First, it reduces the ambient gas densities such that the SN bubble will retain the overpressurized hot gas for a longer time as the radiative cooling rate becomes lower than the energy injection rate of the temporally clustered SNe (Gatto et al. 2015; Kim & Ostriker 2015; Walch et al. 2015; Girichidis et al. 2016; Gentry et al. 2017). Second, it generates coherent gas flows with less momentum cancellation. Therefore, it facilitates the formation of large SN bubbles that can break out of the gaseous disk (Kim et al. 2017; Orr et al. 2022).

The SN clusters are defined by the friends-of-friends method following Smith et al. (2021). This is possible as we record the location and time of each SN event in the simulations. In addition to a linking length of 10 pc, we also impose a linking time of 1 Myr to ensure that the SNe are indeed temporally clustered. The number of SNe in each SN cluster is defined as the clustering number N_cl. We perform the analysis in the "lab" frame and do not correct for the effect of galactic rotation. This means that an SN might drift away from the SN cluster it physically belongs to by more than a linking length due to rotation and thus it would be incorrectly excluded. We argue that this should only be a minor effect as most of our SNe are temporally clustered at much shorter timescales than 1 Myr (see Figure 10).

Figure 7 shows the spatial distributions of the top 10 SN clusters (ranked by N_cl) in different models. Each point represents an SN event, and different colors represent different SN clusters. The average clustering number of the top 10 SN clusters 〈N_cl〉 is shown in each panel.

Zoom In Zoom Out Reset image size

Comparing the fiducial models, gizmo* and arepo* show significantly more clustered SNe than ramses, with 〈N_cl〉 larger by about an order of magnitude. In Lagrangian codes, as gas is able to collapse to densities far above the threshold density, star formation is locally enhanced due to the ${\rho }_{{\rm{g}}}^{1.5}$ dependency. This makes the SNe highly clustered, which collectively generates large SN bubbles and low ${n}_{\mathrm{SN}}$, and leads to bursty star formation. The significantly less clustered SNe in ramses suggests that having a low ${n}_{\mathrm{SN}}$, as we showed in Figure 5, is not a guarantee for efficient feedback—the SNe have to be clustered, too. Meanwhile, the location of SNe in gizmo* and arepo* closely follows the trajectory of galactic rotation, implying that the SN progenitors rarely drift away from their birth clouds and explode in the diffuse medium, which is a potential cause for the low ${n}_{\mathrm{SN}}$. This is perhaps not surprising as we do not resolve the collisional stellar dynamics and neither do we include any subgrid treatment for the so-called "runaway" or "walkaway" stars (de Mink et al. 2014).

To further support our argument that strong SN clustering is indeed caused by locally enhanced SFR, we conduct another numerical experiment, presented in Appendix B, where we force the local SFR to scale linearly rather than superlinearly with gas density. In this case, despite the formation of dense clouds, star formation in the clouds proceeds slowly, resulting in low clustering of SNe and an attendant lack of burstiness (resulting in inefficient outflows—see the following section). In this way, the Lagrangian code behaves similarly to the fiducial Ramses model despite using the fiducial SN delay time of 10 Myr, emphasizing that the root difference is tied ultimately to SN clustering.

On the other hand, the instantaneous SN models (gizmo_tsn0, arepo_tsn0, and ramses_tsn0) represent the extreme case where SN clustering is strongly suppressed. Once star formation occurs, SN feedback immediately stops any further gravitational collapse and the development of any subsequent SNe. This leads to high ${n}_{\mathrm{SN}}$ and thus more rapid energy loss, as well as smaller SN-driven bubbles, resulting in a weaker dynamical impact on the ISM.

With _SF = 100%, SNe are highly clustered in all models. The Lagrangian codes show a similar level of SN clustering compared to their fiducial counterparts (where _SF = 1%). The lack of clusters in the central area of arepo_sfe100 is caused by its strong initial blowout. As opposed to the ramses model, ramses_sfe100 shows significantly more clustered SNe, with 〈N_cl〉 larger by a factor of 30, consistent with its low ${n}_{\mathrm{SN}}$, bursty SFR, and large SN bubbles.

Figure 8 shows N_cl as a function of the median time of the SNe in a cluster (t_cl) in the gizmo* model. Each circle represents an SN cluster. The effect of the SN delay time on clustering is clearly demonstrated: once ${t}_{\mathrm{SN}}$ is switched from 0 to 10 Myr at t = 250 Myr (indicated by the vertical dashed line), SNe immediately become significantly clustered.

Zoom In Zoom Out Reset image size

To demonstrate the development of SN bubbles, Figure 9 shows the SN ambient temperature (${T}_{\mathrm{SN}}$, left panel) and density (right panel) as a function of the normalized time (t − t_cl) for the top four SN clusters (ranked in N_cl) in the gizmo* model. Each point represents an SN event and different colors and symbols represent different clusters. All of these clusters share a similar time evolution: the first SNe occur in cold and dense gas, which gradually heat up and evacuate the gas, eventually saturating at ${T}_{\mathrm{SN}}\sim {10}^{8}$ K and ${n}_{\mathrm{SN}}\sim {10}^{-3}\,{\mathrm{cm}}^{-3}$.

Zoom In Zoom Out Reset image size

To quantify the environmental properties of SN clusters, we calculate the median of the ambient temperature (T_cl) and density (n_cl) of each SN cluster as a function of N_cl in Figure 10. The solid line shows the median values of T_cl or n_cl in each N_cl bin while the shaded area brackets the 25th and 75th percentiles. For gizmo_tsn0 and arepo_tsn0, SNe are forced to occur in cold and dense gas by construction. In contrast, in both gizmo and arepo, as N_cl increases, T_cl increases while n_cl decreases. At N_cl > 20, T_cl becomes higher than 10⁵ K, indicating the development of hot gas in the SN bubbles, which has been shown to be critical for launching galactic outflows (Walch et al. 2015; Girichidis et al. 2016; Hu 2019).

Zoom In Zoom Out Reset image size

Galactic outflows driven by SNe play a critical role in galaxy formation. We characterize the outflows in a similar fashion as in Hu (2019). We define the mass outflow rate as

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{out}}=\displaystyle \sum _{i}\displaystyle \frac{{m}_{{\rm{g}},i}{v}_{{\rm{z}},i}}{{\rm{\Delta }}z}\end{eqnarray} \tag{ 3 }$$

and the energy outflow rate as

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{out}}=\displaystyle \sum _{i}\displaystyle \frac{{m}_{{\rm{g}},i}({v}_{i}^{2}+\gamma {u}_{i}){v}_{{\rm{z}},i}}{{\rm{\Delta }}z},\end{eqnarray} \tag{ 4 }$$

where v_i is the gas velocity, v_z,i is the gas velocity along the vertical direction, γ = 5/3 is the adiabatic index, u_i is the specific internal energy, and Δz = 0.1 kpc is the thickness of the plane where outflows are measured. Here z = 0 corresponds to the midplane of the disk in the initial conditions. The summation is over cells or particles that are (i) located between z = 10 ± 0.5Δz kpc or z = −10±0.5 Δ z kpc and (ii) traveling "outward," i.e., z_i v_z,i > 0.

We define the mass-loading factor,

$$\begin{eqnarray}&&{\eta }_{{\rm{m}}}\equiv \displaystyle \frac{{\dot{M}}_{\mathrm{out}}}{\langle \mathrm{SFR}\rangle },\end{eqnarray} \tag{ 5 }$$

and the energy-loading factor,

$$\begin{eqnarray}&&{\eta }_{{\rm{e}}}\equiv \displaystyle \frac{{M}_{\mathrm{SN}}{\dot{E}}_{\mathrm{out}}}{{E}_{\mathrm{SN}}\langle \mathrm{SFR}\rangle },\end{eqnarray} \tag{ 6 }$$

where ${E}_{\mathrm{SN}}={10}^{51}$ erg is the injection energy per SN and ${M}_{\mathrm{SN}}=100\,{M}_{\odot }$ is the corresponding mass formed in a stellar population per SN. The bracket 〈...〉 represents time-averaging, and we exclude t < 250 Myr to discard the initial transient phase. We use the time-averaged SFR instead of the instantaneous SFR as the normalization factor such that the fluctuations of η_m and η_e are purely due to the fluctuations in the outflows. The time-averaged η_m and η_e, summarized in Table 3, are therefore $\langle {\eta }_{{\rm{m}}}\rangle =\langle {\dot{M}}_{\mathrm{out}}\rangle /\langle \mathrm{SFR}\rangle$ and $\langle {\eta }_{{\rm{e}}}\rangle ={M}_{\mathrm{SN}}\langle {\dot{E}}_{\mathrm{out}}\rangle /({E}_{\mathrm{SN}}\langle \mathrm{SFR}\rangle )$, respectively.

Figure 11 shows ${\dot{M}}_{\mathrm{out}}$ (top panel) and ${\dot{E}}_{\mathrm{out}}$ (bottom panel) as a function of time across the planes of z = ±10 kpc for the Lagrangian models. Instantaneous SN models have more than two orders of magnitude lower outflow rates compared to the fiducial models, demonstrating that SN clustering substantially enhances outflows. Both gadget_tsn0 and gizmo_tsn0 show weak but nonzero outflows at t ∼ 100 Myr as a result of the initial star formation, which is absent in arepo_tsn0. This is likely due to the difference in the initial conditions: arepo_tsn0 includes a low-density background for numerical purposes while gadget_tsn0 and gizmo_tsn0 adopt a vacuum background that presumably facilitates the relatively weak outflows to expand and fill up the vacuum.

Zoom In Zoom Out Reset image size

Interestingly, the outflow rates in the fiducial models show notable differences even though their SFRs, ${n}_{\mathrm{SN}}$, and n_SF are all in good agreement. The difference between arepo* and gadget* is mostly due to their SFRs, which differ by a factor of 3, while their 〈η_m〉 and 〈η_e〉 only differ by a factor of 25% and 80%, respectively. However, gizmo* is a factor of 4–5 lower in 〈η_m〉 and a factor of 10–18 lower in 〈η_e〉 compared to arepo* and gadget*. The origin of this difference is still unclear and requires further investigation in future work.

We now compare outflows between Eulerian and Lagrangian codes in Figure 12. For clarity, we only show the Arepo models as a representation of Lagrangian codes. The ramses model is a factor of 32 lower in 〈η_m〉 and a factor of 100 lower in 〈η_e〉 compared to arepo*, reflecting the fact that SNe are significantly more clustered in arepo*. The outflow rates in Ramses are enhanced with either a larger ${t}_{\mathrm{SN}}$ or a larger _SF, as both would increase the SN clustering. Indeed, ramses_tsn50 and ramses_sfe100 both show comparable outflow rates with arepo*.

Zoom In Zoom Out Reset image size

We conclude that SN clustering plays a fundamental role in driving galactic outflows. Only models with significant SN clustering show 〈η_m〉 larger than unity.

Several recent studies have quantified outflows in terms of η_m and η_e in similar initial conditions of a dwarf galaxy with resolved SN feedback. Our fiducial Lagrangian models (gizmo*, arepo*, and gadget*) are in broad agreement with these studies: Hu (2019) found 〈η_m〉 ∼ 4 and 〈η_e〉 ∼ 0.07 at ∣z∣ = 10 kpc using the Gadget code. Smith et al. (2021) found η_m fluctuating between 1 and 10 and η_e between 0.003 and 0.3 at ∣z∣ = 10 kpc in a comparable model (their PE-SN ) using the Arepo code. Gutcke et al. (2021) found η_m between 5 and 10 and η_e between 10⁻⁴ and 10⁻³ at ∣z∣ = 2 kpc in a comparable model (their Fixed_SN_energy) using the Arepo code. We note that η_m is in better agreement among these studies than η_e, probably because η_e is very sensitive to small differences in the gas velocity (see Equation (4)). Our results confirm these previous findings that the predicted η_m from SN-resolved galaxy-scale simulations are significantly lower than what cosmological simulations commonly assume in their subgrid models. Observations of nearby dwarf galaxies seem to support the low-η_m case (McQuinn et al. 2019).

In a similar setup, Hu et al. (2017) found that a shortened SN delay time (${t}_{\mathrm{SN}}=3$ Myr) has the effect of suppressing the formation of SN bubbles and galactic outflows, which is consistent with our results. However, they did not find an effect of ${t}_{\mathrm{SN}}$ on the burstiness of star formation. Smith et al. (2021) showed that photoionization from massive stars prior to the SN events (the so-called "early feedback") reduces SN clustering and thus significantly suppresses the burstiness of star formation, formation of SN bubbles, and galactic outflows, reducing η_m and η_e by more than an order of magnitude. This is very similar to our instantaneous SN models, which can be viewed as an extreme case of early feedback. Keller & Kruijssen (2022) simulated a Milky Way–like galaxy with m_g ∼ 10⁵ M_⊙, where SN feedback is modeled in a subgrid fashion. Although the simulated galaxy and the numerical treatments are very different, they found that a longer SN delay time enhances galactic outflows by enhancing the clustering of young stars, which is qualitatively consistent with our results.

Kim et al. (2016) have conducted a detailed code comparison study for an isolated Milky Way–like galaxy with significantly more participating codes. They found that different codes generally agree well with each other. In particular, they did not find systematic differences between Lagrangian and Eulerian codes in terms of burstiness and the sizes of SN bubbles. This is probably because they have adopted the simple thermal injection as their subgrid prescription for SN feedback, which is known to be very inefficient at this resolution as most energy would be radiated away without generating much momentum.

As we have shown in Section 3.3, the differences between Lagrangian and Eulerian codes arise near the resolution limit when the Jeans length becomes unresolved. This is perhaps not surprising: Code differences almost by definition have to occur near the resolution limit, as the resolved scales should converge to the physical solutions for any numerically consistent method. However, in our case, the differences at the resolution limit quickly propagate to much larger, well-resolved scales due to clustered star formation and SN feedback.

As we alluded to in Section 3.3, once the gas enters the Jeans length-unresolved regime, both methods can no longer faithfully follow the collapse, and the evolution of gas depends sensitively on numerics such as gravitational softening. That said, it is interesting to ask which method is "less wrong." For a collapsing cloud that becomes unresolved, both the Jeans mass and Jeans length would keep decreasing as the density increases. In other words, the cloud should keep collapsing and fragmenting. Lagrangian codes cannot capture the fragmentation without particle splitting or cell refinement, and they might underestimate the collapsing speed if the adopted softening length is much larger than the cell size. However, the cloud will continue to collapse as expected. In contrast, the collapse in an Eulerian code would halt when the Jeans length becomes unresolved as it requires accreting gas into the minimum cell. In this sense, Lagrangian codes are arguably more accurate than Eulerian codes in the unresolved regime. ¹⁵

The situation becomes more subtle if we consider the physical processes we do not include in this work. In particular, as we do not include feedback from photoionization, we might have overestimated SN clustering and more significantly so in our Lagrangian models. If this is the case, our Eulerian model would more accurately, though coincidentally, predict the SN clustering.

Previous studies have shown that the star formation efficiency only has a weak effect on the galaxy-scale SFR in both Lagrangian and Eulerian codes (e.g., Agertz et al. 2013; Benincasa et al. 2016; Semenov et al. 2017; Hopkins et al. 2018). While this is consistent with our Lagrangian models, it appears to be in conflict with our Eulerian models where increasing _SF leads to a qualitative change in the burstiness of star formation due to enhanced SN clustering. This is probably because SN feedback is unresolved and treated in a subgrid fashion in those studies, which reduces its dynamical impact. Alternatively, as those studies simulated more massive galaxies, the local burstiness might have been averaged out such that the global SFR remains the same. In this case, we expect the outflow rates to increase with _SF.

Instead of assuming a constant _SF, some recent simulations of galaxy formation have adopted a class of subgrid models for star formation that calculates _SF based on local gas properties (Semenov et al. 2016; Kretschmer & Teyssier 2020). These models assume a subgrid log-normal density distribution as predicted by simulations of supersonic turbulence, and models differ in their criteria for the onset of star formation (Krumholz & McKee 2005; Hennebelle & Chabrier 2008; Padoan & Nordlund 2011; Federrath & Klessen 2012; Burkhart 2018; Burkhart & Mocz 2019). Our results suggest that adopting such a subgrid model would have a much stronger effect on Eulerian codes as they are more sensitive to the variation of _SF.

The significant differences between our low- and high-resolution models suggest that we have not yet reached numerical convergence. This is consistent with Smith et al. (2018) whose models with m_g = 20 M_⊙ and m_g = 200 M_⊙ still show significantly different outflow rates. In addition, the convergence study in Hu (2019) with m_g = 1, 5, 25, and 125 M_⊙ showed that convergence was only achieved at m_g = 5 M_⊙ when individual SNe are properly resolved. Therefore, we optimistically expect that our high-resolution models are actually close to convergence.

Many recent cosmological "zoom-in" simulations have adopted an SN feedback model that injects thermal energy for resolved SNe and injects the terminal momentum for unresolved SNe as a subgrid model (e.g., Hopkins et al. 2014, 2018; Agertz & Kravtsov 2015; Kimm et al. 2015; Marinacci et al. 2019). Most of these simulations have resolutions much coarser than our low-resolution models, suggesting that their SN feedback still operates in the subgrid regime where galactic outflows are expected to be underestimated due to the lack of hot gas (Hu 2019). More recent studies such as Agertz et al. (2020), Gutcke et al. (2022), and Calura et al. (2022) have started to resolve individual SNe in cosmological "zoom-in" simulations of dwarf galaxies, which is a promising way forward.

The main caveat of our results is the fact that we do not include pre-SN feedback such as photoionization. As we already discussed, this could lead to overestimated SN clustering, in particular for our Lagrangian models. This is a potential solution for the discrepancies we find between the two types of methods. An interesting follow-up project is therefore to include photoionization and see if the discrepancies can be mitigated. On the other hand, the suppression of SN clustering by photoionization in Smith et al. (2021) could also have been overestimated by the adopted "Strömgren-type" approach where dense clumps are preferentially ionized. Interestingly, Rathjen et al. (2021) conducted ISM-patch simulations using radiative transfer based on an inverse ray-tracing technique and reached the same conclusion that photoionization suppresses SN clustering, which leads to less efficient hot gas generation. This should be investigated by future galaxy-scale simulations using more accurate methods of radiative transfer such as adaptive ray tracing (Emerick et al. 2019). We note that including a solver for radiative transfer does not necessarily resolve the issue, as the important part is to ensure sufficient angular resolution to prevent the dense clumps from being preferentially ionized. This therefore implies a very demanding computational cost.

Another potentially important element is the collisional stellar dynamics. While we adopt a collisionless gravity solver in this work, the stellar dynamics in young star clusters is actually collisional. As our results highlight the importance of SN clustering, it is desirable to include an accurate N-body integrator, such as in Wall et al. (2020) and Rantala et al. (2021), to properly follows the evolution of young star clusters. However, the computational feasibility of incorporating it into galaxy-scale simulations remains to be explored. On the other hand, stellar dynamics can be modeled in a subgrid fashion. Steinwandel et al. (2022) simulated a dwarf galaxy similar to ours, including a subgrid model for runaway massive stars that can enhance 〈η_m〉 by 50% and 〈η_e〉 by a factor of 5.