# Bursty Star Formation in Dwarfs is Sensitive to Numerical Choices in
Supernova Feedback Models

Eric Zhang Department of Physics and Astronomy, University of California,
Riverside, CA, 92507, USA  Laura V. Sales Department of Physics and Astronomy,
University of California, Riverside, CA, 92507, USA  Federico Marinacci
Department of Physics & Astronomy “Augusto Righi”, University of Bologna, via
Gobetti 93/2, 40129 Bologna, Italy  INAF, Astrophysics and Space Science
Observatory Bologna, Via P. Gobetti 93/3, 40129 Bologna, Italy  Paul Torrey
Department of Astronomy, University of Virginia, 530 McCormick Road,
Charlottesville, VA 22903, USA  Mark Vogelsberger Department of Physics, Kavli
Institute for Astrophysics and Space Research, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  Volker Springel Max-Planck-Institut für
Astrophysik, Karl-Schwarzschild-Straße 1, 85740 Garching bei München, Germany
Hui Li  Department of Astronomy, Tsinghua University, Haidian DS 100084,
Beijing, China  Rüdiger Pakmor  Max-Planck-Institut für Astrophysik, Karl-
Schwarzschild-Straße 1, 85740 Garching bei München, Germany  Thales A. Gutcke
Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu,
HI 96822, USA

###### Abstract

Simulations of galaxy formation are mostly unable to resolve the energy-
conserving phase of individual supernova events, having to resort to subgrid
models to distribute the energy and momentum resulting from stellar feedback.
However, the properties of these simulated galaxies, including the morphology,
stellar mass formed and the burstiness of the star formation history, are
highly sensitive to numerical choices adopted in these subgrid models. Using
the SMUGGLE stellar feedback model, we compute idealized simulations of a
$M_{\rm vir}\sim 10^{10}\,\textup{M}_{\odot}$ dwarf galaxy, a regime where
most simulation codes predict significant burstiness in star formation,
resulting in strong gas flows that lead to the formation of dark matter cores.
We find that by varying only the directional distribution of momentum imparted
from supernovae to the surrounding gas, while holding the total momentum per
supernova constant, bursty star formation may be amplified or completely
suppressed, and the total stellar mass formed can vary by as much as a factor
of $\sim 3$. In particular, when momentum is primarily directed perpendicular
to the gas disk, less bursty and lower overall star formation rates result,
yielding less gas turbulence, more disky morphologies and a retention of cuspy
dark matter density profiles. An improved understanding of the non-linear
coupling of stellar feedback into inhomogeneous gaseous media is thus needed
to make robust predictions for stellar morphologies and dark matter core
formation in dwarfs independent of uncertain numerical choices in the baryonic
treatment.

galaxies: structure – galaxies: starburst – galaxies: dwarf – galaxies:
evolution – galaxies: star formation – galaxies: haloes – methods: numerical

Numerics of SNe and Bursty Star Formation Numerics of SNe Feedback

††thanks: NASA Hubble Fellow

##  1 Introduction

Galaxy formation at the scale of dwarfs ($M_{*}\sim
10^{6}-10^{9}\,\textup{M}_{\odot}$) is extremely inefficient, turning only a
small fraction ($\sim 10\%$) of the available gas into stars. One of the most
important processes affecting the formation of stars in dwarf galaxies are
supernovae (SNe), which deposit a significant amount of energy into the
interstellar medium (ISM) (e.g., Larson, 1974; Navarro & White, 1993; Agertz
et al., 2013). This energy deposition is instrumental in regulating the rate
of star formation in galaxies of all masses to the observed levels (e.g.,
Agertz et al., 2010; Guedes et al., 2011; Hopkins et al., 2011; Aumer et al.,
2013). On the scale of dwarfs, supernova energy coupled efficiently to the ISM
may also be responsible for generating fluctuating gravitational potentials
that result in the redistribution of dark matter from the center to the outer
regions of the galaxy and the formation of constant density cores (e.g.,
Navarro et al., 1996a; Pontzen & Governato, 2012; Di Cintio et al., 2014;
Dutton et al., 2016; El-Badry et al., 2016; Read et al., 2016; Tollet et al.,
2016; Freundlich et al., 2020; Jahn et al., 2023; Azartash-Namin et al.,
2024).

In practice, a SNe event can be modeled as a point injection of a large amount
of thermal energy, which expands as a shock into its environment. The first
stage of this expansion is energy-conserving, as the thermal pressure inside
the shock drives the bubble’s expansion according to the Sedov-Taylor solution
(e.g., Taylor, 1950; Sedov, 1959). Once the shock slows down and the driving
pressure becomes comparable to the ambient pressure, the bubble enters the
snow-plow or momentum-conserving phase, where the swept up matter continues to
move without being subject to external forces.

However, it is very difficult for simulations to accurately capture the way in
which SNe explosions expand into the ISM (e.g, Scannapieco et al., 2012).
Firstly, the Sedov-Taylor solution, and related supernova feedback models
(e.g., Chevalier, 1974; Cioffi et al., 1988; Draine & Woods, 1991; Thornton et
al., 1998), invoke a number of simplifying assumptions – in particular, it is
derived based on a homogeneous local ISM at rest. However, neither of these
assumptions are necessarily true for SNe in galaxies; the local ISM may be
inhomogeneous on a wide range of scales, including scales comparable to the
size of the blast wave itself. It will also often not be at rest, and in many
cases we should expect the gas to be moving away from active star-forming
areas, due to previous SNe events or the action of stellar winds. Furthermore,
the resolution of simulations is often insufficient to explicitly capture the
relevant physics of the SNe blastwave. The time resolution needed to resolve
the radiative cooling of very hot gas after a point injection of thermal
energy is difficult to reach in galaxy formation simulations, leading to what
is known as the ‘overcooling problem’ (e.g., Katz, 1992; Smith et al., 2018;
Vogelsberger et al., 2020), in which, due to purely numerical errors, the
thermal energy injected by supernovae dissipates too quickly to expand as a
shock into the ISM. Often, the spatial resolution is also too coarse to
individually resolve the energy-conserving and momentum-conserving phases of
the expanding shock. Typical cooling radii, or size of the bubble at the onset
of the momentum-conserving phase, span from $\sim$ pc-scale to a few hundred
parsecs, depending on the local conditions of the ISM. This is often smaller
than the spatial resolution of galaxy-scale numerical simulations, leaving the
energy-conserving phase unresolved.

As such, SNe feedback is often implemented as a sub-grid process. Methods
include modified thermal injection schemes (e.g., Murante et al., 2010; Dalla
Vecchia & Schaye, 2012), directly injecting momentum into the ISM (e.g.,
Navarro & White, 1993; Mihos & Hernquist, 1994; Vogelsberger et al., 2013;
Hopkins et al., 2018a), temporarily disabling radiative cooling after SNe
events (e.g., Stinson et al., 2006; Agertz et al., 2010; Teyssier et al.,
2013), modeling clusters of multiple SNe collectively as an expanding
superbubble (e.g., Keller et al., 2014; El-Badry et al., 2019), or resorting
to smoothly representing the ISM by means of an effective equation of state
(Springel & Hernquist, 2003). Several current subgrid SNe models based on
direct moment injection have begun to explicitly compensate for resolution
limitations by injecting the terminal momentum of the SNe ejecta after the
unresolved Sedov-Taylor phase, accounting in the calculation for various
dependencies including, local gas density and metallicity (e.g., Kim &
Ostriker, 2015; Pittard, 2019; Hopkins et al., 2018a; Marinacci et al., 2019;
Karpov et al., 2020). Notably, FIRE-3 (Hopkins et al., 2023) adjusts the
momentum injected also for instances when the medium is not at rest, in an
attempt for a more physically-motivated model. In general, higher order
considerations and more detailed sub-grid treatments may be advisable
(Hopkins, 2024).

Recently, some simulations have reached the ultra-high resolution needed to
model SNe via direct injection of thermal energy (e.g., Emerick et al., 2019;
Hu, 2019; Lahén et al., 2020; Gutcke et al., 2021; Deng et al., 2024), and can
thus explicitly resolve the expansion of SNe bubbles into realistic,
inhomogeneous ISM structures. At these resolutions, the momentum injection
into the ISM will naturally arise from the thermal expansion, and will thus
not rely on arbitrary numerical choices. However, these resolutions are not
yet achievable in large-scale simulations. The reality of SNe injection in
galaxy formation simulations is, likely, substantially more complex than
assumed in the prescriptions used, mostly due to the (unresolved)
inhomogeneities in the ISM of galaxies.

For example, even a modestly low density of $n\sim 0.01$ particles per ${\rm
cm}^{3}$ corresponds to a cooling radius of about $200$ pc, which is often
comparable to the size of superbubbles formed by clustered SNe feedback (e.g.,
Li et al., 2024). Studies using small-scale simulations of a more realistic
inhomogeneous ISM (e.g., Martizzi et al., 2015) have shown departures in the
evolution of the SNe remnants from the traditionally assumed sub-grid models
that rest on assumptions of spherically symmetric expansion into an
homogeneous media. The propagation speed and momentum coupling of the SNe
energy vary along low- and high-density channels within the structured ISM,
contributing to the departure from idealized solutions in homogeneous media.
On a macroscopic level, for stratified media (such as the disk of galaxies)
this can result in preferential coupling of the energy/momentum for different
directions, which might not be properly taken into account when enforcing a
symmetric energy/momentum deposition in subgrid models.

It is necessary to produce accurate subgrid models of SNe feedback since many
of the predicted properties of simulated galaxies may show variations
depending on the details of the baryonic treatment and, in particular, the
stellar feedback implementation (Sales et al., 2010; Scannapieco et al.,
2012). This is particularly true in the regime of dwarf galaxies, where the
relatively shallower gravitational potential makes the system more responsive
to the energy and momentum input from feedback. There is consensus in the
field that stellar feedback, and mostly SNe, are capable of regulating star
formation in dwarf galaxies to observed levels, by disrupting cold, dense gas
clumps that would form far too many stars compared to observations if not
regulated by feedback (e.g., see comparison of several codes in Sales et al.,
2022). However, there is less agreement on the level of “burstiness" that is
reasonable to expect in the regime of dwarfs, and whether theoretical models
compare well with observational constraints (e.g., Sparre et al., 2017; Patel
et al., 2018; Emami et al., 2019, 2021).

If star formation is indeed as bursty as predicted by some of the state-of-
the-art numerical simulations, this carries important consequences for the
theoretical predictions on the stellar sizes, morphologies and inner dark
matter distribution in dwarf galaxies. If gas contributes sufficiently to the
local density in the inner regions of halos, which may occur in the case of
gas-rich dwarfs, bursty episodes of star formation can generate substantial
gas flows that create fluctuating gravitational potentials in the inner
regions of low mass halos combined with powerful outflows (e.g., Smith et al.,
2021). These quick (i.e., on sub-dynamical timescales) fluctuations in the
gravitational potential of the galaxy result in a net gaining of energy for
collisionless components like stars and dark matter, making them less tightly
bound, and thus they tend to migrate outwards (Navarro et al., 1996a; Pontzen
& Governato, 2012; Benítez-Llambay et al., 2019). Shallower distributions of
stars and dark matter are thus expected in systems where star formation is
bursty and stellar feedback is efficiently coupled to the surrounding ISM,
including a transformation from cuspy dark matter halos into cored ones.
Through this process, numerical details on the sub-grid modeling of star
formation and stellar feedback have a direct impact on cosmological
predictions and our understanding of CDM as plausible dark matter model.

In this paper, we revisit the assumption of an isotropic deposition of SNe
momentum for sub-grid models of feedback using the multi-phase SMUGGLE model
(e.g., Marinacci et al., 2019) applied in simulations of an isolated dwarf
halo. In Sec. 2 we discuss the simulation set-up and the methods used to
control both: $i)$ directional anisotropy, and $ii)$ cell weighting scheme of
the SNe feedback model. In Sections 3 and 4 we present our main results on how
varying these two assumptions may impact the large-scale dwarf properties,
including their star formation histories, morphologies and inner dark matter
distributions. In Section 5 we summarize our results.

##  2 Description of the Simulations

We use the moving-mesh hydrodynamics code AREPO (Springel, 2010) to run
simulations of an isolated $M_{\rm vir}\sim 10^{10}\,\textup{M}_{\odot}$,
$M_{*}\sim 10^{8}\,\textup{M}_{\odot}$ dwarf galaxy where virial quantities
are defined within the radius that encloses an average density of $200$ times
the critical density of the Universe. Our runs have a reference resolution
with $\sim 6000$ $\textup{M}_{\odot}$ for gas cells and $\sim 25000$
$\textup{M}_{\odot}$ for the dark matter particles. We test convergence of the
results using a $\times 10$ better resolution run in the Appendix.

The initial conditions consist of a gas-rich disk ($M_{\rm gas}\sim 5\times
10^{8}\,\textup{M}_{\odot}$) in hydrostatic equilibrium with the dark matter
halo, with an additional population of old disk stars ($M\sim 7\times
10^{7}\,\textup{M}_{\odot}$) and old bulge stars ($M\sim 5\times
10^{6}\,\textup{M}_{\odot}$). The initial dark matter halo is an NFW (Navarro
et al., 1996b) profile with concentration $c=15$, the initial gas and stellar
disks are in an exponential profile, and the initial stellar bulge is a
Hernquist (1990) profile. This set up is similar to those presented in Hopkins
et al. (2011); Burger et al. (2022a, b); Jahn et al. (2023). The galaxy lies
at the center of a 3D box with 200 kpc on a side. All simulations are run for
at least 2 Gyr, with a snapshot saved every 5 Myr. All galaxy-wide properties
discussed within this work are calculated from all snapshots over the first 2
Gyr, unless stated otherwise.

We briefly summarize the SMUGGLE stellar feedback model as per Marinacci et
al. (2019) (hereafter M+19). It includes various processes of star formation,
stellar feedback and the heating and cooling of gas. Star formation is
primarily modulated by two numerical parameters: the star formation efficiency
$\epsilon_{\rm SF}$, which regulates the speed at which a given gas cells form
stars, and a density threshold $\rho_{\rm th}$, which gas cells are required
to exceed in order to form stars. Chief among the relevant feedback
mechanisms, and the one we will focus on in this paper, is the feedback done
by supernovae.

The main parameters involved in the star-formation and feedback modeling are
summarized in Table 1. Our findings are robust to a range of numerical
choices, which is presented in Appendix A.

Parameter | Value
---|---
Halo Virial Mass | $1\times 10^{10}\,\textup{M}_{\odot}$
Disk Stellar Mass (Initial) | $7\times 10^{7}\,\textup{M}_{\odot}$
Bulge Stellar Mass (Initial) | $5\times 10^{6}\,\textup{M}_{\odot}$
Gas Mass (Initial) | $5\times 10^{8}\,\textup{M}_{\odot}$
Target Gas Mass | 6000 $\textup{M}_{\odot}$
Dark Matter Particle Mass | 25000 $\textup{M}_{\odot}$
SF Efficiency | 0.01
SF Dens. Threshold | 100 $m_{\rm H}/{\rm cm^{3}}$
Table 1: List of the numerical parameters and initial conditions used in the
simulations.

###  2.1 Star Formation

To model star formation, gas cell are converted into star particles,
probabilistically, at a rate consistent with

$\dot{M}_{\rm SF}=\epsilon_{\rm SF}\frac{M_{\rm cell}}{t_{\rm dyn}},$ (1)

where $\epsilon_{\rm SF}$ is the star formation efficiency, which is a tunable
dimensionless factor which we set to 0.01 (M+19, in accordance with Krumholz &
Tan, 2007), $M_{\rm cell}$ is the mass of the gas cell and $t_{\rm dyn}$ is
the dynamical timescale of the gas within the cell, given by:

$t_{\rm dyn}=\sqrt{\frac{3\pi}{32G\rho_{\rm cell}}},$ (2)

where $G$ is the gravitational constant and $\rho_{\rm cell}$ is the density
of the gas cell.

Additionally, gas cells must be sufficiently dense in order to form stars;
this is enforced by allowing only gas cells with a density $\rho_{\rm
cell}>\rho_{\rm th}$ to become a star particle. In the main body of this work,
we set $\rho_{\rm th}$ to be $100$ particles per ${\rm cm}^{3}$. The SMUGGLE
model also requires that the gas is locally gravitationally bound in order to
form stars; this is modeled by a virial parameter criterion for each gas cell,
described in M+19; only gas cells satisfying this criterion are eligible to
become a star particle.

At the time of formation, star particles keep the mass, position, and velocity
of their progenitor gas cells. Each star particle represents a stellar
population, consistent with a Chabrier (2003) initial mass function, which
evolves coevally according to the model described in M+19.111For mass
resolutions within this work, the mass of star particles are sufficiently high
that each particle can be treated as a stellar population that samples well
the IMF.

###  2.2 Supernova Feedback

SNe events occur according to the stellar population evolution model described
in M+19, in which Type II or Ia supernovae occur in a star particle according
to the ages of OB and white dwarf stars within a given stellar population. SNe
are resolved discretely; typically, at the resolution of our simulations, no
more than 1 SN occurs within a star particle at any given timestep.

At each supernova event, a nominal energy of

$E_{\rm tot}=n_{\rm SN}E_{51}\,,$ (3)

is released by a star particle into the neighboring gas particles where
$n_{\rm SN}$ is the number of supernovae occurring within the numerical event,
and $E_{51}=1.0\times 10^{51}$ ${\rm erg}$. This energy is entirely converted
into a scalar momentum, via

$p_{\rm tot}=\sqrt{2m_{\rm tot}E_{\rm tot}},$ (4)

where $m_{\rm tot}$ is the total mass of the supernova ejecta. This linear
momentum is then partitioned among the nearest $N_{\rm ngb}=64$ gas cells,
according to a weight factor $w_{i}$ for each neighboring cell, indexed by
$i$. Then, the (nominal) linear momentum injected into each neighbor gas cell
is

$\tilde{p}_{i}=\frac{p_{\rm tot}w_{i}}{\sum_{i}w_{i}},$ (5)

so that, in principle, the sum of all linear momenta entering the gas cells is
$p_{\rm tot}$. However, a number of numerical ‘post-processing’ procedures are
applied before imparting the momentum to a gas cell. Firstly, due to the
unresolved nature of the energy-conserving phase of the SN blast wave in these
simulations, the linear momentum injected into a gas cell $i$ by a star
particle is boosted222Due to the typical value of $p_{\rm term}$ given in Eq.
7, this value is in practice always greater than 1. (in a manner consistent
with Hopkins et al., 2018a) by a factor

$\beta=\min\left(\sqrt{1+\frac{m_{i}}{\delta m_{i}}},\frac{p_{\rm
term}}{p_{\rm tot}}\right),$ (6)

where $m_{i}$ is the mass of the target gas cell, $\delta m_{i}=m_{\rm
tot}w_{i}/\sum_{i}w_{i}$ is the mass of the ejecta entering the gas cell, and
$p_{\rm term}$ is the terminal momentum of the SNe shock, which is (in
accordance with Cioffi et al., 1988, and used by M+19), given by

$p_{\rm term}=4.8\times 10^{5}\,\textup{M}_{\odot}\,{\rm km}\,{\rm
s}^{-1}\cdot\bar{E}^{13/14}\bar{n}^{-1/7}f(Z)\,,$ (7)

where $\bar{E}$ is the total energy released in units of $10^{51}$ erg,
$\bar{n}$ is the local density estimate in units of particles per ${\rm
cm}^{3}$, and

$f(Z)=\min\left[\left(\frac{Z}{Z_{\odot}}\right)^{-0.14},2\right]^{3/2}\,,$
(8)

where $Z$ is the metallicity. The inclusion of the minimum function in Eq. (6)
ensures that even when the SN momentum is boosted, the total momentum injected
does not exceed the maximum possible momentum injected. The actual linear
momentum injected into a gas cell is then

$p_{i}=\beta\tilde{p}_{i}\,.$ (9)

Furthermore, the nearest $N_{\rm ngb}$ neighboring gas cells are sometimes too
far away from the star particle for momentum injection to be considered
physical. To prevent star particles from being able to affect far-away gas
cells, a ‘superbubble limiter radius’ of $r_{\rm SB}=1024\,{\rm pc}$ is
imposed on all star particles, preventing momentum from being injected into
gas cells further than this radius.333M+19 verifies that varying this radius
within a large range of values does not change the results.

The default choice for the weight is the solid angle of opening $\Omega_{i}$,
defined by

$\Omega_{i}=2\pi\left\\{1-\frac{1}{\left[1+A_{i}/(\pi|{\bf r}_{i}-{\bf
r}_{*}|^{2})\right]^{1/2}}\right\\}\ ,$ (10)

where $A_{i}$ is the cell boundary area. Note that this weighting scheme
follows that in other models, such as FIRE (see Eq. 2 in Hopkins et al.,
2018a). However, other choices for the cell weights are possible, such as the
kernel-weighted cell volume $w_{i}=W_{i}V_{i}$ and the kernel-weighted cell
mass $w_{i}=W_{i}m_{i}$ (where $V_{i}$ is the volume of the target gas cell,
and $W_{i}=W(|{\bf r}_{i}-{\bf r}_{*}|,h)$ is the standard cubic spline SPH
kernel value).

The weight factor $w_{i}$ controls the magnitude of momentum that each gas
cell receives. The direction of this momentum injection can then independently
be determined by the nominal anisotropy factor, $\zeta$. In particular, we
control the vector momentum injected into each gas cell according to an
anisotropic direction vector ${\bf k}$, such that the momentum ${\bf\delta
p_{i}}$ delivered to the gas particle is

${\bf\delta p_{i}}=\frac{p_{i}{\bf k}}{\left|{\bf k}\right|}\ ,$ (11)

where the anisotropic direction vector ${\bf k}$ is defined by

${\bf k}=\left[x_{i}-x_{*}\,,y_{i}-y_{*}\,,\zeta(z_{i}-z_{*})\right]\ ,$ (12)

where $(x_{i},y_{i},z_{i})$ is the position of the target gas particle, and
$(x_{*},y_{*},z_{*})$ is the position of the star particle doing the feedback.
According to this formula, the relative amount of momentum being injected in
the $z$-direction is controlled by the factor $\zeta.$ We show a picture of
this injection scheme in Fig. 1.

Figure 1: Picture of the anisotropic SNe momentum injection scheme, as
parametrized by the nominal anisotropy factor $\zeta$. The star particle
undergoing a SNe event is shown in red, and the surrounding gas cells are
shown as blue circles. The arrows indicate the magnitude and direction of the
momentum injection. A low $\zeta$ value injects momentum primarily within the
disk plane, and a high $\zeta$ injects momentum primarily out of the disk
plane. $\zeta=1$ implies isotropic injection. Note that only the direction of
injected momenta varies with $\zeta$; the magnitude does not.

A higher $\zeta$ means more momentum is being injected in the $z$-direction,
but the normalization factors ensure that the magnitude of momentum released
into a particular gas cell is independent of $\zeta$. $\zeta=1$ implies
isotropic injection; that is, the momentum injected into the cell will always
point radially away from the source star particle. $\zeta$ can only change the
direction of the injected momentum; for any given cell, the magnitude of
momentum injected into it is not affected by the choice of $\zeta$; it is only
controlled by the cell weighting scheme $w_{i}$ and boost factor $\beta$.

As mentioned before, some post-processing calculations may affect the final
momentum injected into the cell after the weight and direction has been
defined. This occurs due to variations in the distribution of gas cells about
star particles, which can affect the boost factor (Eq. 6) and superbubble
limiter. As a result, the nominal anisotropy $\zeta$ may be different from the
effective anisotropy imparted to the individual cells. We quantify the true
anisotropy in momentum injection by an effective anisotropy factor, $\gamma$,
which we define by the ratio of the total momentum injected into gas by SNe
along the $z$-axis (regardless of the positive or negative direction), to the
total momentum injected along the $x$-axis444It would also be valid to use the
$y$ axis for this definition, since it is also within the disk plane; however,
we will choose to use the $x$ axis.,

$\gamma=\left<\frac{\sum_{\rm SNe}\sum_{i}\delta p_{i,z}}{\sum_{\rm
SNe}\sum_{i}\delta p_{i,x}}\right>_{\rm stars}.$ (13)

In practice, this factor does not need to be equal to the nominal anisotropy
factor $\zeta$. However, Fig. 2 indicates that $\gamma$ does in fact scale
well with $\zeta$, so that the degree of anisotropy is still reliably
controlled by $\zeta$. Furthermore, $\gamma$ carries an associated
uncertainty, which is the standard deviation of the quantity averaged over all
star particles in Eq. (13).

Thus, we explore the effect of anisotropy by varying $\zeta$, and we label our
runs by the value of $\zeta$ we employ for that run. However, we will measure
all properties in terms of the effective anisotropy, $\gamma$. Changes in the
cell weighting scheme will also be discussed. Table 2 shows a list of all runs
covered in the main body of this paper, the value of $\zeta$ used, and the
resulting value of $\gamma$ for the corresponding simulation.

Figure 2: The effective anisotropy ($\gamma$, as defined in Eq. 13), and its associated uncertainty, as a function of the nominal anisotropy ($\zeta$, as defined in Eqs. 11, 12). The nominal anisotropy controls the ratio of momentum injected in the polar ($z$) direction to the in-plane ($x,y$) directions, by every star particle, before ‘post-processing’ effects are considered. The effective anisotropy measures the ratio between the total momenta actually injected in these directions. Simulation | Nominal SNe | Effective SNe | SNe Cell Weight | log10 $M_{*}$
---|---|---|---|---
Name | Anisotropy $\zeta$ | Anisotropy $\gamma$ | Scheme $w_{i}$ | Formed [$\textup{M}_{\odot}$]
zeta0.1 (In-Disk Favored) | 0.1 | 0.091 $\pm$ 0.030 | Solid Angle ($\Omega$) | 7.80
zeta0.3 | 0.3 | 0.27 $\pm$ 0.08 | Solid Angle ($\Omega$) | 7.71
zeta1 (Isotropic) | 1.0 | 0.87 $\pm$ 0.20 | Solid Angle ($\Omega$) | 7.71
zeta3 | 3.0 | 2.2 $\pm$ 0.4 | Solid Angle ($\Omega$) | 7.58
zeta10 (Out-of-Disk Favored) | 10.0 | 5.6 $\pm$ 1.1 | Solid Angle ($\Omega$) | 7.38
zeta30 | 30.0 | 13 $\pm$ 3 | Solid Angle ($\Omega$) | 7.32
zeta1-VolWgt | 1.0 | 0.90 $\pm$ 0.23 | Kernel-Weighted Volume | 7.81
zeta10-VolWgt | 10.0 | 5.4 $\pm$ 1.4 | Kernel-Weighted Volume | 7.62
zeta1-MassWgt | 1.0 | 0.95 $\pm$ 0.30 | Kernel-Weighted Mass | 7.36
zeta10-MassWgt | 10.0 | 5.6 $\pm$ 1.9 | Kernel-Weighted Mass | 7.37
Table 2: List of the simulations and the parameters used in each. Simulations
are named by the value of zeta ($\zeta$) and the cell weighting scheme used
for that run, but we will quantify any physical properties by their dependence
on the effective anisotropy $\gamma$, and its associated uncertainty.

In theory, differing distributions of gas cells around star particles
throughout the runs can lead to some variation in the total magnitude of
momentum being injected by a typical star particle between runs. In this case,
we would be unable to conclude that differences in the dwarf properties are
purely due to a change in the directional distribution of momentum. However,
we do find that, for the same choice of cell weighting scheme, $w_{i}$,
varying the anisotropy factor does not change the typical momentum released
and imparted to surrounding cells by a star particle555Since the momentum
injected is measured on a per star-particle basis, it will depend on the
simulation mass resolution as well. We do not consider this dependence in this
work, since we use the same mass resolution for all simulations in the main
body of this work. to any significant degree. However, between different cell
weighting schemes, the total momentum injected can differ by more than half an
order of magnitude. We examine the effect of the cell weighting scheme in more
detail in Section 4.

###  2.3 Other Physical Processes Involved In SMUGGLE

In addition to SNe feedback, the SMUGGLE model includes other forms of stellar
feedback, including OB and AGB stellar winds, radiation pressure, and
photoionization of nearby gas. These modes of stellar feedback can be
important for pre-processing gas before SNe go off, but are energetically
subdominant (M+19, Fig. 13). Also included are radiative heating and cooling
mechanisms of the ISM gas, including cooling of H, He, and metals down to 10
K, self-shielding from ionizing UV radiation, cosmic ray heating, and
photoelectric heating. The parameterization of these additional feedback
processes do not change across different runs, so the only variation between
the runs considered in this work are those of the SNe feedback.

##  3 Anisotropic Momentum Injection

We explore the effects of varying the subgrid SNe momentum injection scheme.
In this section, we focus on the effects of changing the directionality of the
momentum injection. As introduced in Sec. 2 and Eq. (13), the anisotropy in
SNe momentum injection is quantified by $\gamma$, where $\gamma=1$ indicates
isotropic injection, and a $\gamma$ of greater (less) than 1 indicates
momentum primarily being injected perpendicular to (within) the plane of the
disk. We will focus on three runs in this section: zeta0.1 (in-disk favored
injection), zeta1 (isotropic injection, default SMUGGLE choice), and zeta10
(out-of-disk favored injection). Where possible we will refer to specific runs
by their descriptive indicators (i.e., ‘isotropic’ for ‘zeta1’). Each galaxy
is followed for at least 2 Gyr.

###  3.1 Star Formation

Of particular interest is the effect on the star formation (SF) history and
its burstiness. A bursty star formation history can determine the mass, size
and morphology of dwarfs, along with their dark matter distribution (e.g.,
Stinson et al., 2007; Pontzen & Governato, 2012; Di Cintio et al., 2014;
González-Samaniego et al., 2014; El-Badry et al., 2016; Sparre et al., 2017).
We plot the star formation histories of the in-disk favored (zeta0.1),
isotropic (zeta1), and out-of-disk favored (zeta10) injection runs in Fig. 3.

Figure 3: SF histories of the isotropic (purple), in-disk favored (green) and
out-of-disk favored (orange) momentum injection runs. The burstiness of the
SFR depends substantially on the injection model, with the isotropic and in-
disk favored models resulting in more bursty histories than the out-of-disk
favored model.

Qualitatively, in-disk favored (green curve) and isotropic injection (purple
curve) are associated with bursty star formation histories, as seen by the
prevalence of sharp peaks in their SFR curves; contrarily, out-of-disk favored
injection (orange curve) leads to a smoother star formation history.

We quantify the burstiness of a run by the parameter

$B=\frac{\sigma\left(\langle{\rm SFR}\rangle_{\tau_{\rm
avg}}(t)\right)}{\langle{\rm SFR}\rangle},$ (14)

where $B$ stands for the burstiness, $\langle{\rm SFR}\rangle_{\tau_{\rm
avg}}(t)$ denotes the rolling average of the SFR over an interval $\tau_{\rm
avg}$, $\sigma$ denotes the standard deviation of the time series, and
$\langle{\rm SFR}\rangle$ is the average SFR of the simulation. The averaging
over short timescales (5 Myr, 10 Myr, 25 Myr) ensure that rapid random
fluctuations in the SFR, which do not affect large-scale movement of gas due
to feedback, are not counted as additional burstiness.

In Fig. 4, we show the dependence of the burstiness on the anisotropy factor
$\gamma$ for various averaging durations $\tau_{\rm avg}$ = 5 Myr, 10 Myr, and
25 Myr. Using this metric, we find that for any choice of $\tau_{\rm avg}$
within this range, isotropic runs with lower $\gamma$ have bursty SF
histories, and runs with $\gamma>1$ have substantially less bursty SF
histories. These findings are robust to several free parameters of the
simulation including the density threshold for star formation, the star
formation efficiency and the numerical resolution (see Appendix A for
details).

We also show the dependence of the total stellar mass on $\gamma$ in the top
panel of Fig. 4 and find that a preferentially off-plane momentum injection
(higher $\gamma$) is associated with a lower stellar mass. The effect is
noticeable, with the final stellar mass differing by roughly half an order of
magnitude between the highest and lowest anisotropy factors considered. From
this, we conclude that out-of-disk momentum injection is associated with a
less bursty and lower overall rate of star formation.

Figure 4: Top: Total stellar mass formed. Higher-$\gamma$ runs (i.e., runs
with more momentum directed outside the disk), form a lower total stellar
mass. The difference in the stellar mass formed can vary by up to half an
order of magnitude, across all runs considered. Bottom: Burstiness parameter
$B$ as defined by Eq. (14). Higher-$\gamma$ runs tend to be less bursty
according to this parameter, over a wide range of SFR averaging timescales.
The horizontal error bars are the uncertainties in $\gamma$, as discussed in
Eq. (13) and Table 2.

###  3.2 Gas Morphology and Bulk Motions

The increased burstiness seen in runs with isotropic or in-disk favored
momentum injection can be explained by differences in how bulk gas motions are
induced by the SN momentum injection. We consider the bulk gas motion
generated in disk plane and out of the disk plane separately. They are
calculated as

$\dot{M}_{R}=\frac{1}{L}\int{\rm d}m\,{\bf v}\cdot{\bf\hat{R}}\ ,$ (15)

for gas movement within the disk, where $\hat{R}$ refers to the unit vector of
the 2D cylindrical radial coordinate, and

$\dot{M}_{z}=\frac{1}{L}\int{\rm d}m\,{\bf v}\cdot\pm{\bf\hat{z}}\ ,$ (16)

for movement perpendicular to the disk. Here ${\bf v}$ is the velocity of the
gas, the integrals are taken over thin shells of thickness $L=50\,{\rm pc}$,
centered about the surface of a cylinder with diameter $2$ kpc and height
$0.5$ kpc, centered about the center of mass of all star particles. This
cylinder is deeply embedded within the gaseous disk, so that the rates
calculated in Eqs. 15 and 16 trace large-scale flows within the disk, rather
than outflows.

Figure 5: SF history (red) measured within a cylinder of radius of $1$ kpc and
height $\pm 0.25$ kpc from the galactic center, and associated mass flow
history of the isotropic (zeta1, top panel) and out-of-disk favored (zeta10,
bottom panel) runs, measured across the boundary of the same cylinder. Radial
bulk gas motions are shown in blue solid lines, and vertical flows are shown
in green dashed lines. Notably, in the isotropic run, starburst episodes (red
peaks) are preceded by a inward bulk gas motion (blue and green dips) followed
by peaks of outward motion (blue and green peaks) within the disk, consistent
with a more turbulent gas disk. In the in-disk favored run, the vertical gas
cycles are suppressed, and in the out-of-disk favored run, the radial gas
cycles are suppressed. However, only in the out-of-disk favored run is bursty
star formation and the overall mass flows suppressed, indicating that radial
gas flows are more effective at driving turbulence.

Fig. 5 shows the mass flux rate of gas in the plane of the disk (blue) and
along the $z$-direction (green) for the in-disk favored (zeta0.1, top),
isotropic (zeta1, middle) and out-of-disk favored (zeta10, bottom). In the top
and middle panel, we find that in the in-disk favored and isotropic run, a
significant proportion of momentum is injected within the plane of the disk,
leading to large amounts of gas being moved in bulk motions parallel to the
disk plane.

These in-disk bulk gas motions can also drive additional star formation bursts
due to gas compression, leading to additional stars formed in these runs, and
drive instability within the disk, perpetuating the cycle of radially outward
gas motions and bursty star formation. Note the clear correlation between
inward gas motions (dips in the blue curve in Fig. 5), star-formation activity
(peaks in the red curve), and the subsequent outward gas motion (peaks in the
blue curve).

Conversely, in the out-of-disk favored run (bottom panel), the momentum is
primarily injected outside the disk, so that gas compression induced by bulk
radial motions are less prevalent, along with a less bursty star formation
history. This indicates that radial, not vertical, gas flows are primarily
responsible for driving additional star formation and turbulence.
Interestingly, even though the SNe feedback does not lead to bursty SFR or
large-scale gas flows in the out-of-disk favored run, it is still capable of
regulating SFR to similar or lower levels than the isotropic or in-disk
injection runs.

We may also quantify the efficiency of the stars at generating bulk gas flows
via a mass loading factor $\eta$, defined as the ratio of the total gas mass
expelled by SFR-driven flows, to the total stellar mass formed. For a given
run, we calculate

$\eta=\frac{1}{M_{*}}\int\,dt\,\left(\dot{M}_{R,{\rm Out}}+\dot{M}_{z,{\rm
Out}}\right)\,,$ (17)

where $M_{*}$ is the total stellar mass formed within the cylinder, and
$\dot{M}_{R,{\rm Out}}$ and $\dot{M}_{z,{\rm Out}}$ are the quantities in Eqs.
15 and 16 integrated only over outward-moving particles within the shell. The
time integral is done over the total simulation time of 2 Gyr. Calculating
this value666Note that this value should not be compared to typical
measurements of the mass loading factors in the literature since we calculate
$\eta$ for a cylinder deeply embedded in the gas disk for the three runs
pictured in Fig. 5 yields $\eta\sim 40$ for the in-disk favored run, $\eta\sim
60$ for the isotropic run, and $\eta\sim 80$ for the out-of-disk favored run.
This indicates that individual stars in the out-of-disk favored run are more
efficient at pushing and driving gas flows than the in-disk favored and
isotropic runs. Conversely, isotropic or in-disk injection models result in
increased turbulence within the disk which is further exacerbated by the
momentum resulting from a larger number of stars formed than in the out-of-
disk model.

Figure 6: Gas morphology of the in-disk favored (left), isotropic (center),
and out-of-disk favored (right) runs at snapshot time 1.2 Gyr. The box width
of each panel is 8 kpc on a side. Runs with higher injected momentum within
the disk exhibit more diffuse structures and larger gas void regions, but
slightly less turbulence in the vertical directions, than runs with more
injected momentum out of the disk. Figure 7: Stellar morphology of the in-disk
favored (left), isotropic (center), and out-of-disk favored (right) runs for
snapshot at 1.2 Gyr. The stellar half-mass radius about the stellar center of
mass is shown by the white dashed lines. The box width of each panel is 6 kpc
on a side.

The effect of the bulk gas motion cycles is further supported by the resulting
structure of the gas disk, shown in Fig. 6. Significant differences in the
structure can be seen as a result of changing the anisotropy factor. In
particular, in-disk (zeta0.1, left panel) and isotropic (zeta1, middle panel)
momentum injection exhibit a turbulent/disrupted gas morphology, including
large regions or bubbles devoid of gas with radius close to $\sim 1$ kpc in
the most extreme cases. This kind of structure is not present in the case of
out-of-disk favored injection (zeta10, right panel), which shows a smoother
and more regular gas disk without the presence of large gas bubbles.

###  3.3  Stellar Morphology and Kinematics

The effect of different schemes for SNe feedback injection extends also to
noticeable changes in the stellar morphology of the simulated dwarf. Previous
claims in the literature (e.g., Stinson et al., 2007; El-Badry et al., 2016)
have highlighted that bursty star formation histories tend to be associated
with a more diffuse distribution of stars, perhaps also contributing to the
erasing of age and metallicity gradients in simulated dwarf galaxies (e.g.
Mercado et al., 2021)

To quantify the stellar structure in our simulations, we start by showing
projections of the stellar morphology of the in-disk favored (left), isotropic
(middle), and out-of-disk favored (right) runs in Fig. 7. As a crude measure
of how centrally concentrated the stars are in the dwarf, we calculate the
stellar half-mass radius $r_{\rm half}$ as the spherical radius from the
stellar center-of-mass in which half of the stellar mass of the galaxy is
contained. We show $r_{\rm half}$ as a white dotted circle in each of the
panels in Fig. 7. $r_{\rm half}$ does not differ significantly between runs
(we have explicitly checked this as a function of time) and, in general, the
projected stellar structure does not have any obvious indication or trend with
the momentum injection model.

Figure 8: Vertical (green) and in-plane (purple) stellar velocity dispersion
of as a function of anisotropy. Changing the degree of anisotropy does not
affect the vertical velocity dispersion of stars, even when we inject more
momentum in that direction. Conversely, the radial velocity dispersion of
stars is inversely correlated with the anisotropy factor, since in-disk
favored momentum deposition heats the orbits of newly formed stars within the
disk.

However, Fig. 8 shows that the kinematics and orbital properties of the stars
have an imprint of the chosen feedback injection scheme. We plot, for each
simulated galaxy, its velocity dispersion in the 2D radial direction (purple)
and the $z$-direction (green), as a function of its anisotropy factor
$\gamma$, where the mean value and uncertainty are calculated over all
snapshots in the latter 1.5 Gyr of evolution. Preferential injection of
momentum in the plane of the disk (lower-$\gamma$ runs) results in noticeably
higher stellar radial velocity dispersion, presumably as a combination of two
effects: new stars being born from gas with already large radial motions, and
the input of dynamical heating to the population of already formed stars as a
result of radial flows and compression waves.

Interestingly, the fact that the stellar velocity dispersion in the
$z$-direction shows little to no dependency on $\gamma$ (green symbols)
suggests that, perhaps, it is the second of these effects that dominates the
orbital properties of the stars. Indeed, if stars were mainly tracing the
properties of outflowing gas from which they are born, a preferential
injection of momentum out-of-disk (large $\gamma$) should inevitably result in
larger $\sigma_{z}$ for the stars, in clear contrast with our measurements.

###  3.4 Distribution of Dark Matter

The difference in the amount and morphology of the bulk gas motions discussed
in the previous sections has the potential to impact not only the baryonic
morphology of the dwarf, but also its dark matter distribution. Bursty star
formation histories that are associated with locally-dominant gas expulsion
can lead to the transformation of cuspy dark matter halos into lower density
dark matter cores (e.g., Pontzen & Governato, 2012; Benítez-Llambay et al.,
2019). However, models resulting in a smoother star formation history (even
when considering the exact same baryonic physics) could in principle preserve
the characteristic $\rm\Lambda CDM$ cuspy profiles, changing the theoretical
expectations for dwarf galaxies.

Figure 9: Dark matter density profiles between 1 and 2 Gyr. The median density
across snapshots at each radius is plotted as a solid line, and the
interquartile range at each snapshot is shaded. Consistent with the
understanding of core formation, the out-of-disk favored run, which has a non-
bursty SF history, has a dark matter cusp, and the isotropic and in-disk
favored runs, which have bursty SF histories, have dark matter cores. For
comparison, the initial cuspy density profile (measured as the median density
from the stellar center-of-mass over the first 200 Myr, similar across all
runs considered) is shown by the grey dashed line.

To address this, we compare in Fig. 9 the dark matter density profile of our
simulated dwarfs in the in-disk favored (zeta0.1, green), isotropic (zeta1,
red) and out-of-disk favored (zeta10, orange) runs. These profiles have been
calculated as the average dark matter density in the last $1$ Gyr to average
out temporal fluctuations. Thick lines indicate the median while the shaded
regions highlight the interquartile range. We find that lower-$\gamma$ runs,
which are associated with more bursty SFR and in-plane gas flows, tend to have
dark matter cores compared to higher-$\gamma$ runs, for which less bursty SFRs
were measured (see Fig. 4) and that retain a substantially more cuspy profile.
For reference, the dashed gray line indicates the initial profile which
corresponds to a cuspy NFW profile with concentration $c=15$. This correlation
between the mode of star formation and the inner dark matter distribution is
consistent with the expectation from the literature, as discussed above.

Notably, the distinction in the cuspy vs. cored dark matter prediction holds
when considering exactly the same physics of the ISM. While previous work has
already highlighted that some baryonic treatments would result in cores while
others would produce cusps (see e.g., Bose et al., 2019) over similar mass
ranges, the treatment of the baryons from the models considered was different,
opening the possibility that the core formation was purely the result of
including more detailed or more “realistic" star-formation and feedback
models. This degeneracy is underscored by models resolving SN directly but not
producing cores (Gutcke et al., 2022a). In short, our results confirm that not
only the physics and treatment for the ISM play a role, but even the detailed
numerical choices for, e.g., distributing the momentum of the stellar feedback
may cause profound changes in the predictions for the density in the inner
regions of dwarfs.

This also suggests that implementations of “universal" corrections for
modifying the cuspy NFW profile into a lower density core that depend solely
on the stellar mass $M_{*}$ or stellar-to-halo mass ratio $M_{*}/M_{\rm vir}$
(as proposed, e.g., by Di Cintio et al., 2014) might describe the behavior for
a given feedback model or family of baryonic models very well, but are not
necessarily applicable to others. For example, all runs considered in Fig. 9
have a stellar-to-halo mass ratio of $\sim 10^{-2.5}$, a regime where the
inner region is predicted to be cored (with a profile slope of $-0.3$ to
$-0.1$). This is consistent with the isotropic and out-of-disk favored runs at
radii below about $\sim 500$ pc, but the profile slope of the out-of-disk
favored run at this distance is nearly $-1$, consistent with a cusp instead.
We hasten to add, however, that this implication regarding universal
modifications to the dark matter profile still needs further confirmation from
runs within the cosmological set-up.

##  4 Robustness of the results to changes in the cell weighting scheme

As stated in Section 2, it is also possible to change the magnitude of the
momentum distributed across gas cells by changing the weighting function, as
described by the weighting factors $w_{i}$. The default cell weighting scheme
in SMUGGLE is based on the solid angle sub-tended by the gas particles from
the target star, in a way that nearby cells or cells with a large area have
larger weights (see Sec. 2 in Marinacci et al., 2019). This scheme follows
other codes that include a similar weighting, like FIRE-2 (Hopkins et al.,
2018b).

However, other choices of the cell weighting scheme are possible, including
mass-based and volume-based weighting. In a hydrodynamic solver with roughly
equal-mass gas cells, these physically correspond to favoring injection
towards dense gas regions, and towards diffuse gas regions, respectively (see
discussion in Hopkins et al., 2018b). In this section, we explore how these
three choices (solid angle, volume, and mass) of cell weighting scheme
influences the feedback results. The formulae for these weighting schemes are
described in Section 2 (in particular, see Eq. 10 and the subsequent
discussion).

Figure 10: Magnitude of momentum injected by SNe per mass of stars formed, for
each run over 2 Gyr. We calculate the total magnitude in two ways: linear (top
panel, Eq. 18), and in quadrature (bottom panel, Eq. 19). For the same
weighting scheme, there is no dependence of the magnitude of momentum injected
on the anisotropy factor $\gamma$. However, between solid-angle weighting
(red), volume weighting (blue), and mass weighting (green), the magnitude of
momentum injected can differ by as much as an order of magnitude. These
behaviors hold using both calculation methods (although the values do depend
on the calculation method). Figure 11: SFR histories (top) and stellar mass
formed (bottom) for runs with different weighting schemes: solid angle (red),
volume (blue) and mass (green). The left and right panels correspond to
isotropic momentum injection ($\gamma\approx 1$) and out-of-disk favored
momentum injection ($\gamma>1$). In general isotropic injection results on
more bursty SFR histories than out-of-disk favored ones. Solid-angle based
criteria assuming isotropic momentum injection maximizes the burstiness of
star-formation. The total stellar mass formed can vary by factors $2$-$3$
between runs with different numerical choices.

Due to possible differences in the configuration of gas cells around the star
particles, and the boost factor, which is calculated individually for each gas
cell according to their instantaneous density and metallicity (see Eq. 32 in
Marinacci et al., 2019), different weighting schemes will not necessarily
result in the same amount of momentum being injected into the ISM. In fact, we
find that the momentum typically injected by a single star particle does
differ across runs. Because the total momentum injected is a vector quantity,
we define two different ways to quantify it. We calculate it, for each star
particle, using the linear method

$p_{\rm tot,Lin.}=\sum_{\rm SNe}|\delta p_{x}|+\sum_{\rm SNe}|\delta
p_{y}|+\sum_{\rm SNe}|\delta p_{z}|,$ (18)

and the quadrature method

$p_{\rm tot,Quad.}=\sqrt{\left(\Sigma_{\rm SNe}\,|\delta
p_{x}|\right)^{2}+\left(\Sigma_{\rm SNe}\,|\delta
p_{y}|\right)^{2}+\left(\Sigma_{\rm SNe}\,|\delta p_{z}|\right)^{2}},$ (19)

where $p_{x}$, $p_{y}$, and $p_{z}$ are the momenta injected during a given
supernova event in the $x$-, $y$-, and $z$-directions, respectively.

Fig. 10 shows, for each simulation, the mean and standard deviation over all
star particles of the quantity

$\left(\frac{p_{\rm tot}}{M_{\rm init}}\right)_{\rm stars},$ (20)

where $p_{\rm tot}$ is given by either Eq. (18) or (19) depending on the
calculation method, and $M_{\rm init}$ is the initial mass of the star
particle. Note that the calculation reflects the mean and dispersion
normalized by stellar mass formed, so the comparison is meaningful even if the
runs formed different amounts of stellar mass by the end. For our default
choice, opening-angle weight (also referred to as “omega"-weight, red
symbols), we show that the median input as a function of anisotropy (effective
anisotropy, $\gamma$) remains quite constant, suggesting that the typical
momentum input for this weighting scheme is independent of the directionality
($\gamma$) chosen. The trend remains in both, the linear (top panel) or in-
quadrature (bottom) measurement of momentum.

We also compare to a volume weighting scheme (blue) and mass weighting scheme
(green) for the isotropic and out-of-disk momentum injection ($\gamma\sim 1$
and $10$, respectively). We have explicitly checked that the effective
anisotropy of momentum injection is comparable across all weighting schemes
(see Table 2), which then allows for a fair comparison between total momentum
injection at a given $\gamma$. Fig. 10 shows a systematic shift of the
momentum effectively imparted per star depending on the choice of weighting
function, with a volume-weighting criteria returning values slightly below but
comparable to the omega weighting scheme, but a mass weighted injection
results in up to $5$ times less momentum input per star. The effect also seems
independent of $\gamma$, at least for the two explored values here. It is
worth noting, that of all weighting methods, the opening angle criteria
adopted by default in SMUGGLE seems to maximize the momentum injection into
the ISM. Also important is the observation that the choice of cell-weighting
may be well physically motivated, but can result in unwanted pathological
behaviors. For instance, mass-weighting schemes in poorly-resolved disks may
lead to the creation of artificially large bubbles (Torrey et al., 2017),
making such a choice perhaps not preferred for runs with intermediate to low
resolution.

These differences in momentum input have a noticeable effect on the star
formation histories of our simulated dwarfs. We show the evolution of the SF
history (top) and the total stellar mass (bottom) in Fig. 11 for these three
weighting schemes, and choosing isotropic momentum injection (left) and out-
of-disk injection (right). The SMUGGLE default choice, omega-weighted and
isotropic momentum injection, leads to the most bursty star formation of all
runs (red curve, top left panel).

All weighting schemes seem to self-regulate star formation, at least for the
$2$ Gyr explored here, with the volume and mass weighting models predicting a
slight but steady increase of SFR with time. Differences in weighting and
directionality of momentum injection can cause variations by factors $2$-$3$
of the total amount of stars formed, with isotropic injection resulting in
systematically larger stellar masses for all weighting schemes, while at fixed
momentum injection directionality, the mass-weighting scheme results in lower
$M_{*}$.

Interestingly, neither the total rate of star formation nor the stellar mass
are correlated with the magnitude of momentum injected shown in Fig. 10. For
example, the mass-weighted runs have the least momentum injected per star, but
also the lowest rate of star formation, contrary to expectations. Furthermore,
the volume-weighted runs form more stars than the omega-weighted runs, but
have less momentum injected per star particle than the omega-weighted runs.
This shows that the effect of momentum injections on the star formation
history is highly non-linear and complex, i.e., the input of more feedback
momentum does not necessarily translate into less stars being formed.

##  5 Conclusions

Undoubtedly, the assumption of isotropy for the momentum deposition resulting
from a single supernova event in a homogeneous and constant density medium is
well justified. However, the ISM of galaxies is far from this idealized
constant density assumption, and stellar clustering may cause the overlap and
fusion of multiple bubbles leading to a complex system where the assumed
isotropy for energy and momentum deposition no longer holds. The scales over
which these mechanisms occur are currently under-resolved in most galaxy-scale
simulations. Here we explored the systematic effects of lifting the hypothesis
of isotropic bubble expansion using idealized numerical simulations of an
isolated dwarf galaxy run with the SMUGGLE model.

We implemented several different schemes to distribute and couple the momentum
resulting from supernova feedback onto neighboring gas cells. More
specifically, we explored independently the effects of $i)$ changing the
directionality of the momentum injection (lifting the isotropic assumption)
and $ii)$ studying the effect of different weighting schemes (independent of
direction). We label our runs according to the excess of momentum deposited
perpendicular to the disk using the parameter $\gamma$. Isotropic momentum
injection corresponds to $\gamma\approx 1$, while $\gamma>1$ and $\gamma<1$
indicate that the momentum injection is favored perpendicular (vertical) and
parallel (horizontal) to the plane of the disk, respectively.

We can summarize our most important results as follows:

  * •

The degree of anisotropy assumed in SNe momentum injection greatly affects the
burstiness of the star formation history. When momentum is primarily directed
perpendicular to the plane of the disk, bursty star formation is suppressed,
and less total stars are formed, compared to runs in which momentum is
primarily injected isotropically or within the plane of the disk.

  * •

Bursty star formation results from runs where momentum is coupled effectively
within the plane of the disk, generating in-plane bulk gas motion in the
cylindrical radial direction, leading to more gas turbulence within the disk.
The gas structure tends to be more disrupted when momentum is primarily
injected within the disk plane, and less disrupted when momentum is primarily
injected perpendicular to the disk plane.

  * •

The anisotropy of the momentum injection also clearly affects the resulting
stellar morphology of the dwarf. This can be observed in the 2D radial
velocity dispersion of stellar particles being more than double in cases of
in-disk favored and isotropic injection compared to an out-of-disk favored
momentum injection scheme. Stellar orbits become clearly less circular for
models where momentum is primarily coupled in the plane of the disk
($\gamma\leq 1$ models), an effect that we attribute to an increased
gravitational heating from in-plane gas flows in these runs compared to
higher-$\gamma$ models.

  * •

Runs with bursty star formation form dark matter cores, and runs with non-
bursty star formation retain their dark matter cusps. This is consistent with
the existing understanding of the relationship between significant gas flows
(typically associated to bursty star formation) and dark matter core
formation. However, our results also suggest that in "mid-resolution"
simulations that rely on subgrid SNe feedback prescriptions, the formation of
dark matter cores may still depend on numerical choices. Thus, the core
properties such as radial extension and inner density profile are not yet
firmly established as a fundamental prediction of the physics included in a
given simulation.

  * •

Different cell weighting schemes actually lead to different amounts of
momentum being injected, which qualitatively changes the star formation
history and may vary the total amount of stellar mass formed by factors
$2$-$3$ within 2 Gyr of evolution. A weighting scheme based on opening angle,
such as that used in the default version of SMUGGLE, seems to maximize the
momentum input per stellar particle, being $\sim 5$ times larger than, for
example, a scheme based on cell mass alone. As such, the cell weighting scheme
must also be understood as a numerical choice (although it could still be
physically motivated).

Our results imply that numerical choices implicit to the sub-grid modeling of
SNe feedback can play a significant role in the predicted properties of
simulated dwarf galaxies. In particular, for the SMUGGLE code, we demonstrate
that the level of burstiness in the star formation history is particularly
sensitive to the numerical implementation of the underlying feedback physics.
This has consequences for the predicted stellar and gas morphology, with more
bursty runs leading to less rotational support of the gas and stars, along
with a more distorted ISM characterized by large bubbles and in-disk gas
flows. While the specific effects of the weighting scheme will likely vary
from model to model (and depend on details of hydrodynamics solver, numerical
implementation, etc), our study is an important reminder that results from
numerical simulations should be scrutinized not only in light of the physics
included, but also of the numerical choices made, a subject that receives
considerably less attention in the literature.

We stress that, though the parameters studied in this work are numerical, this
does not mean that the numerical choices made have no physical correspondence.
Indeed, the anisotropy parameter and the cell weighting scheme may be
physically motivated choices, that depend on the types of gas that SNe
explosions are expected to couple to most effectively. In reality, the true
nature of how SNe momentum couples to ISM gas is not yet well understood, and
subgrid models must accurately mimic this coupling in order to produce true
predictions from simulations. Thus, in future work it may be instructive to
compare these simulations with models that resolve SN explosions explicitly,
such as LYRA (Gutcke et al., 2021, 2022b) and RIGEL (e.g., Deng et al., 2024).
In such ultra-high resolution simulations, details of the momentum
distribution, including the degree of anisotropy, are expected to be produced
self-consistently by the resolved thermal expansion. Detailed investigations
of the bubble morphology in galaxy simulations with resolved SNe might be able
to better inform subgrid models that rely on numerical parameters. These
parameters may depend on more detailed properties of the ISM such as
directional variations in density, in addition to commonly considered "zeroth-
order" parameters such as the energy released or average local density.

A possible interpretation of the results presented here is that the level of
burstiness in simulated dwarfs could be over-represented at least in some
galaxy formation models. For example, the lack of formation of dense and
compact dwarf galaxies in most codes available today supports the scenario
where more violent star-formation and feedback cycles are present in simulated
dwarfs compared to the real Universe (see discussion in Sales et al., 2022).
Similarly, as reported by El-Badry et al. (2018), simulated dwarfs in FIRE on
the scale of $M_{*}\sim 10^{8}$ \- $10^{9}\;\rm M_{\odot}$ seem less
rotationally supported than suggested by HI observations, also consistent with
the possibility of too bursty star formation rates. In this regard, studies
trying to compare the star formation cycles in observation and simulations may
prove one of the most important ways to inform dwarf galaxy formation models
in simulations (e.g., Patel et al., 2018; Emami et al., 2019, 2021; Mehta et
al., 2023).

Constraining the star formation cycles is especially important not only for
our understanding of dwarf galaxy morphology, but also for the underlying link
between bursty star formation histories, gas flows generated and their impact
on the innermost distribution of dark matter in dwarfs. In our runs with less
bursty SFR, the formation of dark matter cores at the center of dwarfs was
strongly suppressed, meaning that our theoretical expectations of a cusp or a
core at the center of dwarfs depends at least on both, the physics included in
the simulation as well as on the numerical choices made in coarsely
implementing this physics in each code.

## Acknowledgements

The authors would like to thank Phil Hopkins and Julio Navarro for insightful
discussions that helped improved the first version of this draft. EZ and LVS
acknowledge the financial support received through NSF-CAREER-1945310 and NSF-
AST-2107993 grants. FM acknowledges funding by the European Union -
NextGenerationEU, in the framework of the HPC project – “National Centre for
HPC, Big Data and Quantum Computing” (PNRR - M4C2 - I1.4 - CN00000013 – CUP
J33C22001170001). PT acknowledges support from NSF-AST 2346977. HL is
supported by the National Key R&D Program of China No. 2023YFB3002502, the
National Natural Science Foundation of China under No. 12373006, and the China
Manned Space Program through its Space Application System. TAG acknowledges
support by NASA through the NASA Hubble Fellowship grant $\\#$HF2-51480
awarded by the Space Telescope Science Institute, which is operated by the
Association of Universities for Research in Astronomy, Inc., for NASA, under
contract NAS5-26555.

##  Appendix A Robustness of Results to Other Numerical Choices

The changes in the dwarf properties that result from varying the anisotropy of
SNe momentum injection are robust to changes in major simulation parameters,
including the efficiency and density threshold of star formation, and the mass
resolution. To verify this we have run 8 additional simulations of the same
dwarf, in which these simulation parameters were varied. These runs are
summarized in Table 3.

In particular, we find that across a wide range of these parameters, varying
the anisotropy factor $\gamma$ produces the same effect: in-disk favored and
isotropic momentum injection produce a bursty SFR history with periods of
major in-plane gas motion, whereas out-of-disk favored injection produces a
much smoother SFR level with a less turbulent bulk gas flows. We show this
behavior in Fig. 12.

Simulation Name | Nominal SNe Anisotropy $\zeta$ | Max. SFR Efficiency | SFR Threshold [$\rm{m_{H}}/\rm{cm}^{3}$] | Target Gas Mass [$\textup{M}_{\odot}$]
---|---|---|---|---
zeta1-HiEff | 1.0 | 0.9 | 100 | 6000
zeta10-HiEff | 10.0 | 0.9 | 100 | 6000
zeta1-LoThr | 1.0 | 0.01 | 1 | 6000
zeta10-LoThr | 10.0 | 0.01 | 1 | 6000
zeta1-HiEff-LoThr | 1.0 | 0.9 | 1 | 6000
zeta10-HiEff-LoThr | 10.0 | 0.9 | 1 | 6000
zeta1-HighRes | 1.0 | 0.01 | 100 | 500
zeta10-HighRes | 10.0 | 0.01 | 100 | 500
Table 3: List of simulations used in our robustness analysis, and the
parameters used in each. Figure 12: Radial and vertical gas flow rate and
associated SFR for all 8 robustness test runs. In all cases, zeta1 runs have
burstier SFR and more violent gas flow patterns, whereas the zeta10 runs do
not. This is the same behavior as observed in the main body of the paper, so
changing the anisotropy affects the runs in this way regardless of these
numerical choices explored here.

## References

  * Agertz et al. (2013) Agertz, O., Kravtsov, A. V., Leitner, S. N., & Gnedin, N. Y. 2013, ApJ, 770, 25, doi: 10.1088/0004-637X/770/1/25
  * Agertz et al. (2010) Agertz, O., Teyssier, R., & Moore, B. 2010, Monthly Notices of the Royal Astronomical Society, 410, 1391, doi: 10.1111/j.1365-2966.2010.17530.x
  * Aumer et al. (2013) Aumer, M., White, S. D. M., Naab, T., & Scannapieco, C. 2013, Monthly Notices of the Royal Astronomical Society, 434, 3142, doi: 10.1093/mnras/stt1230
  * Azartash-Namin et al. (2024) Azartash-Namin, B., Engelhardt, A., Munshi, F., et al. 2024, arXiv e-prints, arXiv:2401.06041, doi: 10.48550/arXiv.2401.06041
  * Benítez-Llambay et al. (2019) Benítez-Llambay, A., Frenk, C. S., Ludlow, A. D., & Navarro, J. F. 2019, MNRAS, 488, 2387, doi: 10.1093/mnras/stz1890
  * Bose et al. (2019) Bose, S., Frenk, C. S., Jenkins, A., et al. 2019, MNRAS, 486, 4790, doi: 10.1093/mnras/stz1168
  * Burger et al. (2022a) Burger, J. D., Zavala, J., Sales, L. V., et al. 2022a, Phys. Rev. Lett., 129, 191103, doi: 10.1103/PhysRevLett.129.191103
  * Burger et al. (2022b) —. 2022b, MNRAS, 513, 3458, doi: 10.1093/mnras/stac994
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392
  * Chevalier (1974) Chevalier, R. A. 1974, ApJ, 188, 501, doi: 10.1086/152740
  * Cioffi et al. (1988) Cioffi, D. F., McKee, C. F., & Bertschinger, E. 1988, ApJ, 334, 252, doi: 10.1086/166834
  * Dalla Vecchia & Schaye (2012) Dalla Vecchia, C., & Schaye, J. 2012, MNRAS, 426, 140, doi: 10.1111/j.1365-2966.2012.21704.x
  * Deng et al. (2024) Deng, Y., Li, H., Liu, B., et al. 2024, arXiv e-prints, arXiv:2405.08869, doi: 10.48550/arXiv.2405.08869
  * Di Cintio et al. (2014) Di Cintio, A., Brook, C. B., Macciò, A. V., et al. 2014, MNRAS, 437, 415, doi: 10.1093/mnras/stt1891
  * Draine & Woods (1991) Draine, B. T., & Woods, D. T. 1991, ApJ, 383, 621, doi: 10.1086/170820
  * Dutton et al. (2016) Dutton, A. A., Macciò, A. V., Dekel, A., et al. 2016, MNRAS, 461, 2658, doi: 10.1093/mnras/stw1537
  * El-Badry et al. (2019) El-Badry, K., Ostriker, E. C., Kim, C.-G., Quataert, E., & Weisz, D. R. 2019, MNRAS, 490, 1961, doi: 10.1093/mnras/stz2773
  * El-Badry et al. (2016) El-Badry, K., Wetzel, A., Geha, M., et al. 2016, ApJ, 820, 131, doi: 10.3847/0004-637X/820/2/131
  * El-Badry et al. (2018) El-Badry, K., Bradford, J., Quataert, E., et al. 2018, MNRAS, 477, 1536, doi: 10.1093/mnras/sty730
  * Emami et al. (2019) Emami, N., Siana, B., Weisz, D. R., et al. 2019, ApJ, 881, 71, doi: 10.3847/1538-4357/ab211a
  * Emami et al. (2021) Emami, N., Siana, B., El-Badry, K., et al. 2021, ApJ, 922, 217, doi: 10.3847/1538-4357/ac1f8d
  * Emerick et al. (2019) Emerick, A., Bryan, G. L., & Mac Low, M.-M. 2019, MNRAS, 482, 1304, doi: 10.1093/mnras/sty2689
  * Freundlich et al. (2020) Freundlich, J., Dekel, A., Jiang, F., et al. 2020, MNRAS, 491, 4523, doi: 10.1093/mnras/stz3306
  * González-Samaniego et al. (2014) González-Samaniego, A., Colín, P., Avila-Reese, V., Rodríguez-Puebla, A., & Valenzuela, O. 2014, ApJ, 785, 58, doi: 10.1088/0004-637X/785/1/58
  * Guedes et al. (2011) Guedes, J., Callegari, S., Madau, P., & Mayer, L. 2011, ApJ, 742, 76, doi: 10.1088/0004-637X/742/2/76
  * Gutcke et al. (2021) Gutcke, T. A., Pakmor, R., Naab, T., & Springel, V. 2021, MNRAS, 501, 5597, doi: 10.1093/mnras/staa3875
  * Gutcke et al. (2022a) —. 2022a, MNRAS, 513, 1372, doi: 10.1093/mnras/stac867
  * Gutcke et al. (2022b) Gutcke, T. A., Pfrommer, C., Bryan, G. L., et al. 2022b, ApJ, 941, 120, doi: 10.3847/1538-4357/aca1b4
  * Hernquist (1990) Hernquist, L. 1990, ApJ, 356, 359, doi: 10.1086/168845
  * Hopkins (2024) Hopkins, P. F. 2024, arXiv e-prints, arXiv:2404.16987, doi: 10.48550/arXiv.2404.16987
  * Hopkins et al. (2011) Hopkins, P. F., Quataert, E., & Murray, N. 2011, MNRAS, 417, 950, doi: 10.1111/j.1365-2966.2011.19306.x
  * Hopkins et al. (2018a) Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018a, MNRAS, 477, 1578, doi: 10.1093/mnras/sty674
  * Hopkins et al. (2018b) —. 2018b, MNRAS, 480, 800, doi: 10.1093/mnras/sty1690
  * Hopkins et al. (2023) Hopkins, P. F., Wetzel, A., Wheeler, C., et al. 2023, MNRAS, 519, 3154, doi: 10.1093/mnras/stac3489
  * Hu (2019) Hu, C.-Y. 2019, MNRAS, 483, 3363, doi: 10.1093/mnras/sty3252
  * Jahn et al. (2023) Jahn, E. D., Sales, L. V., Marinacci, F., et al. 2023, MNRAS, 520, 461, doi: 10.1093/mnras/stad109
  * Karpov et al. (2020) Karpov, P. I., Martizzi, D., Macias, P., et al. 2020, ApJ, 896, 66, doi: 10.3847/1538-4357/ab8f23
  * Katz (1992) Katz, N. 1992, ApJ, 391, 502, doi: 10.1086/171366
  * Keller et al. (2014) Keller, B. W., Wadsley, J., Benincasa, S. M., & Couchman, H. M. P. 2014, MNRAS, 442, 3013, doi: 10.1093/mnras/stu1058
  * Kim & Ostriker (2015) Kim, C.-G., & Ostriker, E. C. 2015, ApJ, 802, 99, doi: 10.1088/0004-637X/802/2/99
  * Krumholz & Tan (2007) Krumholz, M. R., & Tan, J. C. 2007, ApJ, 654, 304, doi: 10.1086/509101
  * Lahén et al. (2020) Lahén, N., Naab, T., Johansson, P. H., et al. 2020, ApJ, 891, 2, doi: 10.3847/1538-4357/ab7190
  * Larson (1974) Larson, R. B. 1974, MNRAS, 169, 229, doi: 10.1093/mnras/169.2.229
  * Li et al. (2024) Li, C., Li, H., Cui, W., et al. 2024, MNRAS, doi: 10.1093/mnras/stae797
  * Marinacci et al. (2019) Marinacci, F., Sales, L. V., Vogelsberger, M., Torrey, P., & Springel, V. 2019, Monthly Notices of the Royal Astronomical Society, 489, 4233, doi: 10.1093/mnras/stz2391
  * Martizzi et al. (2015) Martizzi, D., Faucher-Giguère, C.-A., & Quataert, E. 2015, Monthly Notices of the Royal Astronomical Society, 450, 504, doi: 10.1093/mnras/stv562
  * Mehta et al. (2023) Mehta, V., Teplitz, H. I., Scarlata, C., et al. 2023, ApJ, 952, 133, doi: 10.3847/1538-4357/acd9cf
  * Mercado et al. (2021) Mercado, F. J., Bullock, J. S., Boylan-Kolchin, M., et al. 2021, MNRAS, 501, 5121, doi: 10.1093/mnras/staa3958
  * Mihos & Hernquist (1994) Mihos, J. C., & Hernquist, L. 1994, ApJ, 425, L13, doi: 10.1086/187299
  * Murante et al. (2010) Murante, G., Monaco, P., Giovalli, M., Borgani, S., & Diaferio, A. 2010, Monthly Notices of the Royal Astronomical Society, 405, 1491, doi: 10.1111/j.1365-2966.2010.16567.x
  * Navarro et al. (1996a) Navarro, J. F., Eke, V. R., & Frenk, C. S. 1996a, MNRAS, 283, L72, doi: 10.1093/mnras/283.3.L72
  * Navarro et al. (1996b) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996b, ApJ, 462, 563, doi: 10.1086/177173
  * Navarro & White (1993) Navarro, J. F., & White, S. D. M. 1993, MNRAS, 265, 271, doi: 10.1093/mnras/265.2.271
  * Patel et al. (2018) Patel, S. G., Kelson, D. D., Diao, N., Tonnesen, S., & Abramson, L. E. 2018, ApJ, 866, L21, doi: 10.3847/2041-8213/aae524
  * Pittard (2019) Pittard, J. M. 2019, MNRAS, 488, 3376, doi: 10.1093/mnras/stz1885
  * Pontzen & Governato (2012) Pontzen, A., & Governato, F. 2012, MNRAS, 421, 3464, doi: 10.1111/j.1365-2966.2012.20571.x
  * Read et al. (2016) Read, J. I., Agertz, O., & Collins, M. L. M. 2016, MNRAS, 459, 2573, doi: 10.1093/mnras/stw713
  * Sales et al. (2010) Sales, L. V., Navarro, J. F., Schaye, J., et al. 2010, MNRAS, 409, 1541, doi: 10.1111/j.1365-2966.2010.17391.x
  * Sales et al. (2022) Sales, L. V., Wetzel, A., & Fattahi, A. 2022, Nature Astronomy, 6, 897, doi: 10.1038/s41550-022-01689-w
  * Scannapieco et al. (2012) Scannapieco, C., Wadepuhl, M., Parry, O. H., et al. 2012, MNRAS, 423, 1726, doi: 10.1111/j.1365-2966.2012.20993.x
  * Sedov (1959) Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics
  * Smith et al. (2021) Smith, M. C., Bryan, G. L., Somerville, R. S., et al. 2021, MNRAS, 506, 3882, doi: 10.1093/mnras/stab1896
  * Smith et al. (2018) Smith, M. C., Sijacki, D., & Shen, S. 2018, MNRAS, 478, 302, doi: 10.1093/mnras/sty994
  * Sparre et al. (2017) Sparre, M., Hayward, C. C., Feldmann, R., et al. 2017, MNRAS, 466, 88, doi: 10.1093/mnras/stw3011
  * Springel (2010) Springel, V. 2010, MNRAS, 401, 791, doi: 10.1111/j.1365-2966.2009.15715.x
  * Springel & Hernquist (2003) Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289, doi: 10.1046/j.1365-8711.2003.06206.x
  * Stinson et al. (2006) Stinson, G., Seth, A., Katz, N., et al. 2006, MNRAS, 373, 1074, doi: 10.1111/j.1365-2966.2006.11097.x
  * Stinson et al. (2007) Stinson, G. S., Dalcanton, J. J., Quinn, T., Kaufmann, T., & Wadsley, J. 2007, ApJ, 667, 170, doi: 10.1086/520504
  * Taylor (1950) Taylor, G. I. 1950, Proc. R. Soc. Lond. A, 201, 159
  * Teyssier et al. (2013) Teyssier, R., Pontzen, A., Dubois, Y., & Read, J. I. 2013, MNRAS, 429, 3068, doi: 10.1093/mnras/sts563
  * Thornton et al. (1998) Thornton, K., Gaudlitz, M., Janka, H. T., & Steinmetz, M. 1998, ApJ, 500, 95, doi: 10.1086/305704
  * Tollet et al. (2016) Tollet, E., Macciò, A. V., Dutton, A. A., et al. 2016, MNRAS, 456, 3542, doi: 10.1093/mnras/stv2856
  * Torrey et al. (2017) Torrey, P., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2017, MNRAS, 467, 2301, doi: 10.1093/mnras/stx254
  * Vogelsberger et al. (2013) Vogelsberger, M., Genel, S., Sijacki, D., et al. 2013, MNRAS, 436, 3031, doi: 10.1093/mnras/stt1789
  * Vogelsberger et al. (2020) Vogelsberger, M., Marinacci, F., Torrey, P., & Puchwein, E. 2020, Nature Reviews Physics, 2, 42, doi: 10.1038/s42254-019-0127-2