# On the impact of runaway stars on dwarf galaxies with resolved interstellar
medium

Ulrich P. Steinwandel,1 Greg L. Bryan,2,1 Rachel S. Somerville,1 Christopher
C. Hayward1 and Blakesley Burkhart1,3
1Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue,
New York, NY 10010
2Department of Astronomy, Columbia University, 550 West 120th Street, New
York, NY 10027, USA
3Department of Physics and Astronomy, Rutgers University, 136 Frelinghuysen
Road, Piscataway, NJ 08854, USA E-mail<EMAIL_ADDRESS>(UPS)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

About ten to 20 percent of massive stars may be kicked out of their natal
clusters before exploding as supernovae. These “runaway stars” might play a
crucial role in driving galactic outflows and enriching the circumgalactic
medium with metals. To study this effect, we carry out high resolution dwarf
galaxy simulations that include velocity kicks to massive stars above 8 M⊙. We
consider two scenarios, one that adopts a power law velocity distribution for
kick velocities, resulting in more stars with high velocity kicks, and a more
moderate scenario with a Maxwellian velocity distribution. We explicitly
resolve the multi-phase interstellar medium (ISM), and include non-equilibrium
cooling and chemistry. We adopt a resolved feedback scheme (Griffin) where we
sample individual massive stars from an IMF and follow their radiation input
as well as their the SN-feedback (core-collapse) channel at the end of their
lifetime. In the simulations with runaway stars, we add additional (natal)
velocity kicks that mimic two and three body interactions that cannot be fully
resolved in our simulations. We find that the inclusion of runaway or walkaway
star scenarios has an impact on mass, metal, momentum and energy outflows as
well as the respective loading factors. We find an increase in mass, metal and
momentum loading by a factor of 2-3, whereas we find an increase in the mean
energy loading by a factor of 5 in the runaway case and a factor of 3 in the
walkaway case. We conclude that the inclusion of runaway stars could have a
significant impact on the global outflow properties of dwarf galaxies.

###### keywords:

methods: numerical – galaxies: dwarf – galaxies: formation – galaxies:
evolution – galaxies: ISM – ISM: jets & outflows

††pubyear: 2022††pagerange: On the impact of runaway stars on dwarf galaxies
with resolved interstellar medium–On the impact of runaway stars on dwarf
galaxies with resolved interstellar medium

##  1 Introduction

The $\Lambda$CDM paradigm has been extremely successful at predicting and
reproducing observations of structure formation on large scales (greater than
a few Mpc). However, since the first attempts to model galaxy formation within
this paradigm, it has been clear that _feedback_ from massive stars,
supernovae, and active black holes is critical for reproducing many observed
properties of galaxies within this paradigm (e.g Somerville & Davé, 2015; Naab
& Ostriker, 2017, for reviews). There has been great progress recently in
developing large scale simulations that reproduce many key observations of
galaxy populations (e.g. Hirschmann et al., 2014; Schaye et al., 2015;
Vogelsberger et al., 2014; Pillepich et al., 2018; Nelson et al., 2018). In
addition, in the past decade, numerical simulations finally succeeded in
producing thin, rotation dominated disks with rather inefficient baryon-to-
star conversion (e.g. Guedes et al., 2011; Agertz et al., 2013; Marinacci et
al., 2014; Hopkins et al., 2014, 2018).

However, all of these simulations to one degree or another achieve this
success by implementing somewhat ad hoc or phenomenological “sub-grid”
recipes, as they are not able to explicitly resolve many of the physical
processes that actually drive feedback. For example, large-volume simulations
typically have a spatial resolution of a few hundred parsecs or larger, and
baryonic mass resolution of $10^{5}\mbox{${\rm M}_{\odot}$}$ or greater (see
figure 1 of Nelson et al., 2019 for a recent summmary). Thus these simulations
typically do not resolve the multi-phase interstellar medium (ISM), nor
blastwaves from (individual) supernova explosions, among many other relevant
processes. One of the biggest open questions in galaxy formation theory
currently is how “small scale” processes (that occur below the explicit
resolution of the simulation) may impact larger scale and _global_ properties
of galaxies.

A new generation of simulations are achieving parsec to sub-parsec spatial
resolution and mass resolution comparable to that of single stars. These
simulations can explicitly resolve the multiphase ISM, the Sedov-Taylor phase
of individual supernova explosions, and the multiphase character of outflows,
and are therefore able to implement the physics of star formation and stellar
feedback in a more explicit manner. Several groups have carried out such
simulations for sub-galactic regions (Kim & Ostriker, 2015; Kim et al., 2017;
Kim et al., 2020; Walch et al., 2015), typically representative of Solar
Neighborhood or Milky Way conditions. Recently, a few groups have been able to
achieve this “individual star” resolution and implement these more explicit
feedback approaches within simulations of _entire dwarf galaxies_ with
accurate non-equilibrium heating and cooling (Hu et al., 2017; Emerick et al.,
2019; Lahén et al., 2019b; Smith et al., 2021), which is a significant step
forward. In this work, we leverage this relatively new capability to study a
process that is unresolved in galaxy scale simulations, but may play an
important role in shaping “macroscopic” galaxy properties: runaway stars.

Massive stars are predominantly born in clusters. Simulations of individual
star clusters (e.g. Oh & Kroupa, 2016) have shown that some fraction of
massive stars can receive kicks large enough to move them out of their natal
clouds before exploding as supernovae. These kicks can be imparted by
gravitational interactions within a binary star system (sometimes involving a
third star), or when one of the stars in a binary explodes as a SNae. In both
cases, one of the massive stars can be accelerated above the escape velocity
of the binary system (or even the star cluster itself). These unbound stars
can easily travel at velocities $>30$ km s-1 and can thus travel the typical
length scale of a GMC of a few $10$ parsec on a timescale of $~{}1-2$ Myr.
However, we note that it is not clear yet which of the scenarios for the
formation of runaway stars is favoured (Renzo et al., 2019; Fujii & Portegies
Zwart, 2011).

It is extremely challenging to include these processes in a galaxy scale
simulation, for several reasons. First, one would have to model the binary
stellar intial mass function (IMF), which has significant uncertainties (e.g.
Kroupa, 1995; Malkov & Zinnecker, 2001). Second, and perhaps more challenging,
one would have to model two and three body interactions of single stars.
Galaxy formation simulations are not well equipped to carry out orbital
integration of multiple stellar systems, due to the nature of the gravity
calculation which is often performed with a tree or particle mesh. Hence, the
gravity solvers of many astrophysical simulation codes operate in the limit of
a collisionless system, while this is no longer a good description of few body
interactions in star clusters. At the cloud scale, there are attempts to
circumvent this weakness (e.g. Hirai et al., 2021; Fujii et al., 2021; Grudić
et al., 2022; Guszejnov et al., 2022). However, it has also been pointed out
recently that there might be an issue with the underlying Schmidt-type star
formation recipes often applied in galaxy scale simulations which makes it
intrinsically hard to disperse star clusters on the correct timescales (e.g.
Smith et al., 2021; Hislop et al., 2022), a problem which would likely become
worse when orbits are integrated using direct summation algorithms.

Why might we expect that runaway stars could have a significant effect on
galaxy scale properties? Supernova (SNe) feedback is believed to have the
strongest impact on providing the necessary mid-plane pressure to drive a
galactic outflow. (Martizzi et al., 2015; Kim & Ostriker, 2015; Walch et al.,
2015; Gatto et al., 2017; Naab & Ostriker, 2017; Steinwandel et al., 2020).
However, recently several groups have put forward the idea of galactic winds
that are driven by cosmic-rays (e.g. Buck et al., 2020; Girichidis et al.,
2016; Girichidis et al., 2022; Hopkins et al., 2021a; Hopkins et al., 2021b;
Hopkins et al., 2021c, d, e). For SN-feedback, the crucial condition to drive
outflows is the formation of a hot volume filling ISM phase (e.g. Naab &
Ostriker, 2017; Steinwandel et al., 2020). An important factor in this picture
is the host environment of the SNe (e.g. Hu et al., 2016, 2017, 2019; Lahén et
al., 2019a; Lahén et al., 2019b; Gutcke et al., 2021; Smith et al., 2021;
Hislop et al., 2022). SNe that explode in the cold neutral medium (CNM) or the
cold molecular medium are inefficient in establishing the hot phase of the ISM
as cooling times are short (a few 100 years) and the remnants remain small
($\sim 10$ pc) while in the energy conserving Sedov-Taylor phase. This makes
it relatively unlikely that subsequent supernovae will explode in the bubble
of previous SNe, so that no efficient outflows can be launched from the cold
ISM due to the formation of superbubbles (e.g. Fielding et al., 2017; Orr et
al., 2021, 2022). However, the situation is different in the diffuse ISM at
lower densities and higher temperatures, where SNe remnants occupy a larger
volume fraction. Cooling times are much longer, and a hot phase is naturally
established, due to the larger occupation fraction of the remnants which makes
superbubble formation and subsequent breakout much easier.
There are a limited number of physical processes that can lead to lower
ambient media preceding a SN-event, of which some are driven by the host star
itself (e.g. Forbes et al., 2016; Haid et al., 2016; Hu et al., 2017; Haid et
al., 2018; Emerick et al., 2019; Kim et al., 2019; Smith et al., 2021;
Lancaster et al., 2021a, b, c, for an overview of main and post main sequence
stellar winds, photo-ionisation, photoelectric heating). As discussed above,
runaway stars are another mechanism that could potentially allow massive stars
to travel out of their dense natal gas clouds, such that the SNe potentially
explode in a warmer, lower density environment, enhancing their ability to
drive large scale outflows (e.g. Naab & Ostriker, 2017).

The impact of runaway stars has recently been studied in conditions similar to
the Milky Way in Kim & Ostriker (2018) and Andersson et al. (2020). However,
there has been no detailed study of this issue in the context of dwarf
galaxies, which is the goal of this paper. In this paper, we create 4
simulations of isolated dwarf galaxies similar to the observed Wolff-Lundmark-
Moyette (WLM) system at the outer edge of the Local Group (baryon mass of
$2\cdot 10^{7}$ M⊙ and dark matter mass of $2\cdot 10^{10}$ M⊙). We use the
Griffin(Galaxy Realizations with the feedback from individual massive stars)
model for single star formation and resolved feedback (e.g. Lahén et al.,
2019a; Lahén et al., 2019b; Steinwandel et al., 2020; Hislop et al., 2022). We
include a sub-grid treatment for runaway stars based on existing results for
the runaway population from simulations and observations (e.g. Eldridge et
al., 2011; Oh & Kroupa, 2016). We note that our approach includes only the
natal kicks, neglecting the gravitational interaction of the star with the
environment over its lifetime, and therefore may represent an upper limit on
the effects of including such stars in galaxy simulations.
The remainder of the paper is structured as follows: In Section 2 we introduce
the code, present the Griffin model for single star formation and resolved
feedback, and explain in detail how we incorporate runaway stars into our
model. In Section 3 we discuss the results with a specific focus on how
runaway stars may alter outflow properties of dwarf galaxy systems and impact
the ISM and CGM-ISM interface. In Section 4 we discuss our results and compare
to previous simulations in the literature. In Section 5 we summarise our
results.

##  2 Numerical Methods

In this section we discuss the details of our numerical modelling and present
the initial conditions of our simulations. The results are part of the Griffin
project 111https://wwwmpa.mpa-garching.mpg.de/~naab/griffin-project/.

###  2.1 Simulation Code

For all simulations presented we use the Tree-SPH code P-Gadget3 (Springel,
2005; Hu et al., 2014) to carry out the gravity and hydrodynamics
calculations. However, we made significant improvements by updating the
baseline SPH solver to incorporate the meshless finite mass (MFM) method for
solving the fluid flow following the methods discussed in Hopkins (2015) and
Gaburov & Nitadori (2011), using the MFM implementation of Steinwandel et al.
(2020). MFM represents a fully consistent second order method in space and
time and we compute fluid fluxes by solving the one-dimensional Riemann-
problem between single fluid tracers by adopting an HLLC-Riemann solver that
accounts for the appropriate reconstruction of the contact wave. The face area
is defined over the smoothing-length weighted quadrature point between the
single fluid tracers (Steinwandel in prep.).
The code is then coupled to the cooling and feedback network, utilised in the
Griffin-framework (e.g. Hu et al., 2016, 2017, 2019; Lahén et al., 2019b;
Lahén et al., 2019a; Steinwandel et al., 2020; Hislop et al., 2022).

###  2.2 Cooling, heating and star formation

To model the multiphase ISM, we self-consistently model the non-equilibrium
formation and destruction of molecular hydrogen and CO, following the cooling
and chemistry networks developed in Glover & Mac Low (2007) and Glover & Clark
(2012), where we also account for the self-shielding of molecular gas based on
a column density computation on a healpix sphere. In earlier versions of our
code this was achieved by the TreeCol method. In the latest version this is
done with a healpix expansion using the c library chealpix. The details of the
cooling and chemistry network can be found in Hu et al. (2017) and Steinwandel
et al. (2020).
Star formation is modelled by utilising an IMF-sampling approach where we
sample all the stars more massive than 4 M⊙ (base resolution of the
simulation) as single stars. The details of the method can be found in Hu et
al. (e.g. 2017); we have adopted minor updates here that will be presented in
greater detail in Steinwandel (in prep.).

###  2.3 Stellar feedback

We follow several feedback channels for the massive stars that form within the
simulation, including the photoelectric heating, the UV-radiation and the
photoionising radiation, as well as the metal enrichment from AGB and subgrid
stellar winds. However, the core of the stellar feedback scheme is the
resolved SN-feedback mechanism developed by Hu et al. (2017) and Steinwandel
et al. (2020), that allows for a self consistent build-up of the hot phase and
the momentum of individual SN-events in a resolved Sedov-Taylor phase (see Hu
et al., 2021; Hu et al., 2022). Hence, all the outflow properties in these
simulations are self-consistently driven from the hot phase of the ISM without
the need for any “ad-hoc” tuning of mass loading factors as adopted in lower
resolution galaxy formation simulations. Therefore, our simulations allow for
a detailed study of the origin and driving of large scale galactic winds.

We include the UV radiation of single stars following the implementation of Hu
et al. (2017). This implementation makes the assumption of a (locally)
optically thin medium, which is a good approximation for dwarf galaxies.
Hence, we can simply sum over the radiation field for each star using an
inverse square law. The lifetime functions and effective surface temperatures
are taken from Georgy et al. (2013). The UV luminosity is obtained from the
stellar library BaSel (Lejeune et al., 1997, 1998; Westera et al., 2002) in
the energy band $6-13.6$ eV, which dominates the photo-electric heating rate
in the ISM. We note that one star particle can contain more then one star in
this respect as we trace the UV-radiation from all stars with masses larger
than $2$ M⊙. However, for star particles more massive than $4$ M⊙ the stars
fully represent single stars. We integrate the radiation contributions while
walking the gravity tree. Hence, this is effectively a ray-tracing approach
modulo the angular resolution of the tree. This can be straightforwardly
improved in future implementations by either including higher order moments in
the tree or modifying the tree’s opening angle (i.e. running a separate
gravity tree for the gravity force and UV-radiation computation) (e.g. Grond
et al., 2019, for one approach to improve our method).
We follow photoelectric heating rates from Bakes & Tielens (1994); Wolfire et
al. (2003) and Bergin et al. (2004) at the rate:

$\displaystyle\Gamma_{\mathrm{pe}}=1.3\cdot 10^{-24}\epsilon
DG_{\mathrm{eff}}n\ \ \mathrm{erg}\ \mathrm{s}^{-1}\mathrm{cm}^{-3},$ (1)

with the hydrogen number density $n$, the effective attenuation radiation
field $G_{\mathrm{eff}}=G_{0}\exp{-1.33\cdot 10^{-21}DN_{\mathrm{H,tot}}}$,
the dust-to-gas ratio D, the total hydrogen column density
$N_{\mathrm{H,tot}}$ and the photoelectric heating rate $\epsilon$ that can be
determined via:

$\displaystyle\epsilon=\frac{0.049}{1+0.004\psi^{0.73}}+\frac{0.037(T/10000)^{0.7}}{1+2\cdot
10^{-4}\psi}$ (2)

where $\psi$ is denoted as $\psi=G_{\mathrm{eff}}\sqrt{T}/n^{-}$, with the
electron number density $n^{-}$ and the gas temperature $T$.

We include a treatment for photoionising radiation since previous work (e.g.
Hu et al., 2017; Emerick et al., 2019; Steinwandel et al., 2020; Smith et al.,
2021), has revealed that photoionisation feedback can have a significant
impact on ISM properties and star formation rates. Ideally, we would utilize a
proper radiative transfer scheme for this purpose, but this is prohibitively
computationally expensive. Thus we model the photoionising feedback via a
Stroemgren approximation locally around each photoionising source in the
simulation, similar to the approach used in Hopkins et al. (2012) and later
improved by Hu et al. (2017).
Every star particle that fulfils the criterion of having one or more
constituent stars with mass m${}_{\mathrm{IMF}}>8$ M⊙ is treated as a
photoionising source in our simulations. Thus all of the neighbouring gas
particles are labelled as photoionised within the radius RS around the star
and all of the surrounding molecular and neutral hydrogen is destroyed and
transferred to H+ (we set the ionising fraction to 0.9995 in our chemical
model). Alongside this procedure we heat the surrounding gas within Rs to 104
K, which is the temperature in HII-regions, by providing the necessary energy
input from the ionising source. This procedure is straightforward if RS is
known (i.e. as long as the gas density is more less constant). However, in
practice this will rarely be the case, and we obtain $R_{S}$ iteratively,
where we use the classic Stroemgren-radius as an “educated guess”.
We include the feedback from Type II (CCSN) supernovae in our simulations. The
supernova feedback routines presented in this section are the core feedback
mechanism in our current galaxy formation and evolution framework. Every
supernova distributes $10^{51}$ erg into the nearest 32 neighbours. We follow
metal enrichment using the yields of Chieffi & Limongi (2004), who present
stellar yields for stars from zero metallicity to solar metallicity for
progenitors ranging from 13 M⊙ to 35 M⊙, from which we interpolate our metal
enrichment routines. When a star explodes in a CCSN the mass is added to the
surrounding $32$ (one kernel) gas particles with the respective metal mixture
that is expected from the above yields. It is important to include the correct
metal yields in order to obtain the correct cooling properties of the ambient
medium and thus in the end the correct density and phase structure of the ISM.
However, as the amount of metal mass that is distributed into the neighbours
of the star particle can be much larger than the mass of one gas particle, we
split the gas particles after the enrichment has taken place if the mass
exceeds the initial particle mass resolution of the simulation by a factor of
two, to avoid numerical artefacts. The split particles inherit all the
physical properties that were present in the parent particle including
specific internal energy, kinetic energy, metallicity and chemical abundances.
The particle mass is split by a factor of two and the position is offset by
one fifth of the parent particle smoothing length in a random (isotropic)
direction to avoid the overlap of the two newly spawned particles. The
spawning process requires a new domain decomposition and a re-build of the
gravity tree which makes the code somewhat slower.

Table 1: Overview of the physics variations adopted in our different simulations. For convenience we follow similar naming conventions as those adopted in Hu et al. (2017). Name | Core Collapse Supernovae | Photoionising radiation | Photoelectric heating | runaway stars
---|---|---|---|---
WLM-fiducial | ✓ | ✓ | ✓ | X
WLM-RunP | ✓ | ✓ | ✓ | power law distribution
WLM-RunM | ✓ | ✓ | ✓ | Maxwellian distribution
WLM-inplane | ✓ | ✓ | ✓ | maxwellian (sampling only for vx and vy)

###  2.4 Treatment of Runaway stars

While the spatial and mass resolution of our simulation is relatively high by
current standards, we are not yet capable of accurate modeling of two, three
and few-body interactions needed to self-consistently account for the velocity
offset from which runaway O/B stars originate. We note that attempts can be
made to self-consistently model these phenomena in the near future in galaxy
scale simulations employing the methods presented in Rantala et al. (2017) and
Rantala et al. (2020) (higher order symplectic), and there are recent
demonstrations of similar methods on the cloud scale for higher order hermite
integration of star clusters (e.g Grudić et al., 2022; Guszejnov et al.,
2022). However, currently it is not possible to accurately compute the forces
for two or three body encounters in a galaxy scale simulation that also takes
a detailed ISM model into account. A first order approximation for the effect
of runaway stars can be obtained by constructing a “subgrid” model for the
velocity kicks for the massive O/B stars that become runaway stars. In the
following we present two models, one motivated by simulations for a strong
runaway case and another one which represents a weaker walkway case.

####  2.4.1 The Runaway Case: Power law velocity distribution

In the strong runaway model we assume that the velocity distribution follows a
simple power law

$\displaystyle f_{\mathrm{vel}}\propto v^{-\beta},$ (3)

with the velocity power law slope $\beta$, which is dependent on several
factors. One popular scenario for runaway star formation occurs when a binary
system undergoes gravitational scattering processes with a third star. Thus it
is straightforward to assume that $\beta$ is mainly driven by the number of
binary systems in the star cluster and the relaxation time of the star cluster
itself. Moreover, in reality, star clusters can be disrupted (e.g. due to
internal supernova events or tidal interactions with other star clusters)
which is also important for setting $\beta$. All these processes can lead to
velocity kicks following the distribution of equation 3. However, there would
be a time delay that comes from the direct N-body interactions which is not
captured by our approach. This is a clear limitation of applying random
(isotropic) velocity kicks. To accurately model star cluster formation and
evolution it is necessary to solve the direct N-body system, which would give
a reliable estimate for $\beta$ (Eldridge et al., 2011; Perets & Šubr, 2012;
Oh & Kroupa, 2016; Renzo et al., 2019). In detail, the velocity distribution
may depend on the mass of the star as well as the properties of the star
cluster. Alternatively, one could follow an empirically motivated approach and
adopt an estimate for $\beta$ directly from observations (e.g. Blaauw, 1961;
Comerón & Pasquali, 2007; Gies, 1987; Hoogerwerf et al., 2000, 2001). In our
first runaway model we will follow a theoretically motivated approach by
adopting estimates from Oh & Kroupa (2016). This will allow for a direct
comparison to Andersson et al. (2020), who investigated the effects of runaway
stars in galaxy simulations of the Milky Way.
We start our modelling by building a subgrid model for runaway stars in which
we adopt a power law velocity distribution following eq. 3 with a slope
$\beta=1.8$. To do so we slightly modify our IMF-sampling approach and add
velocity kicks to stars that are more massive than $8M_{\odot}$. First, we
have to sample the velocity distribution of equation 3. This can be
accomplished by drawing a random number $p$ and obtaining the absolute value
for the velocity kick of particle $i$ in the following manner:

$\displaystyle v_{i}=v_{\mathrm{min}}\cdot(1-p)^{-\frac{1}{\beta-1}}.$ (4)

We adopt a minimum value $v_{\mathrm{min}}=3$ km s-1 for direct comparison
with Andersson et al. (2020). We truncate the distribution at $385$ km s-1.
However, this truncation process introduces a complication. If we draw a
random number between zero and one and we calculate $v_{i}$ following equation
4, then it can happen that we naturally sample values above $385$ km s-1. We
want to truncate the sampling because this will give as a distribution where
14 percent of the stars obtain a velocity kick of more than $30$ km s-1, which
is in good agreement with observations. The sampling of values above $385$ km
s-1 can be related to a certain upper limit pmax for the random number $p$,
which is given by

$\displaystyle
p_{\mathrm{max}}=1-\left(\frac{3}{385}\right)^{\frac{4}{5}}\approx 0.9744.$
(5)

We thus only accept the random number $p$ if it is smaller or equal than
$p_{\mathrm{max}}$. This ensures that we obtain a velocity power law
distribution with a slope of $\beta=1.8$ and a $14$ percent fraction of stars
with velocity kicks of more than $30$ km s-1.

####  2.4.2 The Walkaway Case: Maxwellian velocity distribution

The second scenario that we adopt for velocity kicks is designed to mimic the
more conservative case of so called walkaway stars, where the bulk of the
stars in the velocity distribution is moving at a speed lower than 10 km s-1,
following a velocity distribution that is consistent with a Maxwell-Boltzmann
distribution:

$\displaystyle p(v)\propto
v^{2}\exp{\left(-\frac{v^{2}}{2\sigma^{2}}\right)}.$ (6)

where $\sigma$ is the width of the distribution. This is motivated based on
previous work, both observational and theoretical, that shows some indication
for a Maxwellian distribution (e.g. Eldridge et al., 2011; Silva & Napiwotzki,
2011). However, we note that we take a rather conservative approach that is
more consistent with previous research such as Renzo et al. (2019) which is
suggesting that the bulk of runaway stars can actually be classified as
walkaway stars. While such a distribution is a little harder to sample than
the power law from the previous section, it becomes almost trivial when
considering that the distribution of equation 6 is the three dimensional
representation of the Maxwell-Boltzmann distribution integrated over solid
angle in spherical coordinates. Hence, this distribution in Cartesian
coordinates consists of three (independent) Gaussian distributions. The latter
can be straightforwardly sampled by a Gaussian random number generator, which
is available in the gsl package as the Ziggurat-algorithm, and is the fastest
way of sampling the distribution of equation 6. We note that we explored
different options to sample this distribution such as the Box-Muller method
that can obtain Gaussian random numbers by utilising a linear combination of
four uniform random numbers. In practice the difference is marginal, but the
Ziggurat-algorithm is much faster. We note that for this model we also adopt a
probability for a star to get any velocity kick. This is motivated by the fact
a fraction of stars will be born as single stars and a fraction will be born
as binary stars. In this model we assume that only binary stars can receive a
velocity kick and we adopt a binary star fraction of 50 percent222The choice
of 50 per cent is rather ad-hoc since most the exact fraction of
runway/walkaway stars is not exact.. We note that we did not assume this in
the first model, to have a model close to the results of Oh & Kroupa (2016)
and the implementation of Andersson et al. (2020). Note that we do not
explicitly sample the binary IMF, which is a clear limitation of this model.
Since this model is supposed to represent the more conservative walkaway case
we adopt a relatively low value of $\sigma=32$ km s-1.

###  2.5 Initial conditions and simulations

In our set of simulations we use the initial conditions of an isolated dwarf
galaxy with a baryon mass of $2\cdot 10^{7}$ M⊙ and a dark matter mass of
$2\cdot 10^{10}$ M⊙, which is an analogue of a classic Wolff-Lundmark-Moyette
(WLM) system. The mass resolution is set to $4$ M⊙. The force softening for
this system is $0.1$ pc. While our choice for the gravitational softening is
rather aggressive we note that extensive testing of different values of this
parameter has little effect on global galactic properties such as morphology,
star formation, and outflow rates. However, we note that a smaller choice of
the softening-length than $0.1$ leads to sub-optimal gravity time stepping
behaviour due to runaway star formation in dense regions that results in the
gravitational ejection of single stars to arbitrarily large velocity. This
effect is mainly significant in runs without PI-radiation (not presented
here). The softening is roughly comparable to the size of the kernel that we
observe in the dense regions of the ISM during the simulation. The initial
conditions are generated with the code presented in Springel et al. (2005)
often referred to as “MakeDiskGal” or “Makegalaxy”.
In total we carry out four simulations WLM-fiducial, WLM-RunP (power law
runaway stars), WLM-RunM (Maxwellian runaways), and WLM-inplane (runaways only
in xy plane) summarised in Table 1. The first case WLM-fiducial represents a
baseline case without runaway stars. The second case WLM-RunP represents a
strong runaway case. The third case WLM-RunM represents a weaker runaway
scenario, in which stars do not receive such high velocity kicks as in the
WLM-RunP scenario. The final case WLM-inplane samples the same power-law
velocity distribution as WLM-RunP, but kicks are only imparted in the x-y
plane. This allows us to test the impact of runaways migrating within the
plane of the galaxy versus escaping above or below the disk.

###  2.6 Definitions for outflow rates, inflow rates, loading factors and
metal enrichment factors

We follow the definitions for outflow rate, inflow rate, mass loading and
enrichment factors of previous work, specifically, the definitions presented
in Hu et al. (2019). Generally, galactic wind quantities can be split into
mass-flux (m), momentum-flux (p) and energy-flux (e) defined via:

$\displaystyle{\bf\mathcal{F}}_{\mathrm{m}}=\rho\textbf{v},$ (7)
$\displaystyle{\bf\mathcal{F}}_{\mathrm{p}}=\rho\textbf{v}\cdot\textbf{v}+P,$
(8) $\displaystyle{\bf\mathcal{F}}_{\mathrm{e}}=(\rho
e_{\mathrm{tot}}+P)\textbf{v},$ (9)

with $\rho$, $v$ and $P$ given as the fundamental fluid variables density,
velocity and pressure. At all times $e_{\mathrm{tot}}$ is defined as the the
total energy per unit mass $e_{\mathrm{tot}}=0.5\cdot\mathbf{v}^{2}+u$, where
$u$ is the internal energy per unit mass. In standard code units $u$ takes the
units $kms^{-2}$. We note this to avoid confusion with converting these units
to proper cgs units in the sections below.

In this work, if not stated otherwise, all the flow rates are measured at a
height of $1$ kpc above the midplane of the galaxies in a plane-parallel patch
with thickness $dr$, where $dr=0.1$ kpc.

We will adopt the following definitions for in and outflows rates, with the
total flow rates defined via
$\dot{M}=\dot{M}_{\mathrm{out}}-\dot{M}_{\mathrm{in}}$,
$\dot{p}=\dot{p}_{\mathrm{out}}-\dot{p}_{\mathrm{in}}$ and
$\dot{E}=\dot{E}_{\mathrm{out}}-\dot{E}_{\mathrm{in}}$. Outflow is defined by
a positive value of the radial velocity given via
$\mathbf{v\cdot\hat{n}}=v_{r}>0$, and inflow is defined as the material in the
patch with negative radial velocity $v_{r}<0$. Hence, we can write down the
discrete inflow and outflow rates for mass, momentum and energy:

$\displaystyle\dot{M}_{\mathrm{out}}=\sum_{i,v_{i,r}>0}\frac{m_{i}v_{i,r}}{dr},$
(10)
$\displaystyle\dot{M}_{\mathrm{in}}=-\sum_{i,v_{i,r}<0}\frac{m_{i}v_{i,r}}{dr},$
(11)
$\displaystyle\dot{p}_{\mathrm{out}}=\sum_{i,v_{i,r}>0}\frac{m_{i}[v_{i,r}^{2}+(\gamma-1)u_{i}]}{dr},$
(12)
$\displaystyle\dot{p}_{\mathrm{in}}=-\sum_{i,v_{i,r}<0}\frac{m_{i}[v_{i,r}^{2}+(\gamma-1)u_{i}]}{dr},$
(13)
$\displaystyle\dot{E}_{\mathrm{out}}=\sum_{i,v_{i,r}>0}\frac{m_{i}[v_{i,r}^{2}+\gamma
u_{i}]v_{i,r}}{dr},$ (14)
$\displaystyle\dot{E}_{\mathrm{in}}=-\sum_{i,v_{i,r}<0}\frac{m_{i}[v_{i,r}^{2}+\gamma
u_{i}]v_{i,r}}{dr},$ (15)

where we use $P=\rho u(\gamma-1)$ as an equation of state. Furthermore, we
will compute metal outflow rates via:

$\displaystyle\dot{M}_{\mathrm{out}}=\sum_{i,v_{i,r}>0}\frac{Z_{i}m_{i}v_{i,r}}{dr},$
(16)
$\displaystyle\dot{M}_{\mathrm{in}}=-\sum_{i,v_{i,r}<0}\frac{Z_{i}m_{i}v_{i,r}}{dr},$
(17)

where $Z_{i}$ is the metallicity of particle $i$. It is useful to investigate
these parameters normalised to a set of reference quantities, defined as the
so-called loading factors for mass ($\eta_{\mathrm{m}}$), metals
($\eta_{\mathrm{Z}}$), momentum ($\eta_{\mathrm{p}}$) and energy
($\eta_{\mathrm{e}}$) and given by:

  1. 1.

outflow mass loading factor:
$\eta_{\mathrm{m}}^{\mathrm{out}}=\dot{M}_{\mathrm{out}}/\overline{\mathrm{SFR}}$,

  2. 2.

inflow mass loading factor:
$\eta_{\mathrm{m}}^{\mathrm{in}}=\dot{M}_{\mathrm{in}}/\overline{\mathrm{SFR}}$,

  3. 3.

momentum outflow loading factor:
$\eta_{\mathrm{p}}^{\mathrm{out}}=\dot{p}_{\mathrm{out}}/(p_{\mathrm{SN}}R_{\mathrm{SN}})$,

  4. 4.

momentum inflow loading factor:
$\eta_{\mathrm{p}}^{\mathrm{in}}=\dot{p}_{\mathrm{in}}/(p_{\mathrm{SN}}R_{\mathrm{SN}})$,

  5. 5.

energy outflow loading factor:
$\eta_{\mathrm{e}}^{\mathrm{out}}=\dot{E}_{\mathrm{out}}/(E_{\mathrm{SN}}R_{\mathrm{SN}})$,

  6. 6.

energy inflow loading factor:
$\eta_{\mathrm{e}}^{\mathrm{in}}=\dot{E}_{\mathrm{in}}/(E_{\mathrm{SN}}R_{\mathrm{SN}})$,

  7. 7.

outflow metal loading factor:
$\eta_{\mathrm{Z}}^{\mathrm{out}}=\dot{M}_{\mathrm{Z,out}}/(m_{\mathrm{Z}}R_{\mathrm{SN}})$,

  8. 8.

inflow metal loading factor:
$\eta_{\mathrm{Z}}^{\mathrm{in}}=\dot{M}_{\mathrm{Z,in}}/(m_{\mathrm{Z}}R_{\mathrm{SN}})$,

where we adopt $E_{\mathrm{SN}}=10^{51}$ erg, p${}_{SN}=3\cdot 10^{4}$ M⊙ km
s-1, the supernova rate R${}_{SN}=\overline{\mathrm{SFR}}/(100M_{\odot})$ and
the metal mass (IMF-averaged) $m_{\mathrm{Z}}=2.5M_{\odot}$ which we take for
best comparison with previous dwarf galaxy studies from Hu et al. (2019). We
compute the mean star formation rate $\overline{\mathrm{SFR}}$ over the time
span of one Gyr, which is the total run time of the simulation. Additionally,
we define metal-enrichment factors that not only quantify the amount of metals
transported in the outflows, but can quantify to what degree metals are over-
or under-abundant in galaxy outflows relative to the ISM. This can be achieved
by normalising the metal outflow rate to the mass outflow rate, weighted by
the background metallicity of the galactic ISM,
Z${}_{\mathrm{gal}}=0.1Z_{\odot}$:

  1. 1.

outflow enrichment factor:
$y_{\mathrm{Z}}^{\mathrm{out}}=\dot{M}_{\mathrm{Z,out}}/(\dot{M}_{\mathrm{out}}Z_{\mathrm{gal}})$,

  2. 2.

inflow enrichment factor:
$y_{\mathrm{Z}}^{\mathrm{in}}=\dot{M}_{\mathrm{Z,out}}/(\dot{M}_{\mathrm{in}}Z_{\mathrm{gal}})$.

##  3 Results

In this section we will study the outflow properties of all of our dwarf
galaxy simulations with a specific focus on the time evolution of the outflow
rates and wind loading factors.

Figure 1: Visualisation of the the face-on stellar surface density (left), the
gas surface density (center) and the pressure (right), for our isolated dwarf
galaxy simulations, for the runs WLM-fiducial (first row), WLM-RunP (second
row), WLM-RunM (third row) and WLM-inplane (fourth row). We find only very
minor differences between the different runs in terms of the morphology of the
stellar structure and the ISM.

Figure 2: We visualise the edge-on stellar surface density (left), the gas
surface density (center) and the temperature (right), for our isolated dwarf
galaxy simulations, for the runs WLM-fiducial (first row), WLM-RunP (second
row), WLM-RunM (third row) and WLM-inplane (fourth row). We find only very
minor differences between the different runs in terms of the morphology of the
stellar structure and the ISM. However, the stellar disc seems to be slightly
“puffed up” in the runs WLM-RunP and WLM-RunM compared to the runs WLM-
fiducial and WLM-inplane.

Figure 3: We show the star formation rate for the runs WLM-fiducial (black),
WLM-RunP (blue), WLM-RunM (red) and WLM-inplane (orange) in the top panel, as
well as the total stellar mass of the systems in the bottom panel as a
function of time. We find a slight boost of the star formation rate in the
runaway star models. This may be due to the reduction of feedback from PI-
radiation as some massive stars quickly leave their natal clouds. However, the
time averaged star formation rate is only weakly influenced by this effect due
to the self-regulating nature of the star formation and feedback cycle. Figure
4: Distribution of the density of the gas in which SNe explode, for all runs
WLM-fiducial (black), WLM-RunP (blue), WLM-RunM (red) and WLM-inplane
(orange). For all runs we took a statistically large enough sample of 3000 SNe
explosions to allow for direct comparison. The data are binned in 128 equally
spaced logarithmic bins in the range from $10^{-4}$ to $10^{1}$ cm-3. Runaway
stars have only a weak effect on the distribution function but in the run WLM-
RunP there is a slight excess of stars exploding in lower density
environments. However, this is within the model scatter, which implies that
major changes in global properties are driven by the altitude above or below
the disc at which the SNe explode and not the redistribution of the SNe in the
galactic plane.

Figure 5: We show the time averaged mass weighted density-temperature phase
space for the runs WLM-fiducial (top), WLM-RunP (centre) and WLM-RunM
(bottom). All plots are averaged over a time scale of 500 Myr. The runs WLM-
RunP and WLM-RunM show excess mass in the hot phase of the ISM. However, since
the hot phase typically does not contain most of the mass, the trend is more
apparent in the volume weighted version of this figure (Fig. 6).

Figure 6: We show the time averaged volume weighted density-temperature phase
space for the runs WLM-fiducial (top), WLM-RunP (centre) and WLM-RunM
(bottom). All plots are averaged over a time scale of 500 Myr. It is apparent
that the inclusion of runaway stars can have a significant effect on
establishing the hot, volume filling phase of the ISM. Figure 7: We show the
volume filling factor of the hot gas (above a temperature of $3\times 10^{4}$
K) as a function of time for all four models WLM-fiducial (black), WLM-RunP
(blue), WLM-RunM (red), and WLM-inplane (orange). One can clearly see that the
runaway models show an increase in the volume filling factor of the hot phase.
As pointed out in previous studies of Hu et al. (2017) and Steinwandel et al.
(2020), the volume filling factor is critical for driving outflows in
galaxies.

Figure 8: We show the outflow velocity as a function of the sound speed in
logarithmic space for the models WLM-fiducial (top) and WLM-RunP. The colour
code indicates the averaged outflow mass flux rate over a time span of 500
Myr. While these indicate that the bulk of the mass in the outflow is
transported by the cold phase, we find a significant increase in mass
transport in the hot phase of the wind, that extends to higher outflow
velocities and sound speed for the simulation WLM-RunP (bottom).

Figure 9: We show the outflow velocity as a function of the sound speed in
logarithmic space for the models WLM-fiducial (top) and WLM-RunP. The colour
code indicates the averaged outflow energy flux rate over a time span of 500
Myr. The energy flux carried by the phase at high velocities and high sound
speeds is significantly increased when runaway stars are included in the run
WLM-RunP (bottom).

Figure 10: We show the mass outflow rate (top left), the metal outflow rate
(top right), the momentum outflow rate (bottom left) and the energy outflow
rate (bottom right) as a function of time for all four simulated models, where
black indicates the fiducial model WLM-fiducial, blue represents the power law
model WLM-RunP and red the Maxwellian runaway model WLM-RunM. The orange line
shows the runaway model WLM-inplane in which we apply the velocity kicks only
along the vx and the vy directions. The impact of runaway stars can be clearly
seen in a boost of the mass, the metal and the momentum outflow rate by a
factor of around two (comparing blue lines to black lines). In the case of the
energy outflow rate the boost is even more significant, around a factor seven
to eight.

Figure 11: We show the mass outflow loading (top left), the metal outflow
loading (top right), the momentum outflow loading (bottom left) and the energy
outflow loading (bottom right) as a function of time for all four simulated
models, where black indicates the fiducial model WLM-fiducial, blue represents
the power law model, WLM-RunP and red the Maxwellian runaway model WLM-RunM.
The orange line shows the runaway model WLM-inplane in which we apply the
velocity kick only along the vx and the vy directions. While the trend of
increased mass and metal outflow rates is reduced when expressed in terms of
the normalized loading factors, the trend of increased momentum and energy
loading remains roughly the same as that deduced from the results of Fig. 10.

Figure 12: We show the mass inflow rate (top left), the metal inflow rate (top
right), the momentum inflow rate (bottom left) and the energy inflow rate
(bottom right) as a function of time for all four simulations WLM-fiducial
(black), WLM-RunP (blue), WLM-RunM (red) and WLM-inplane. We find increased
inflow rates for all the models that include runaway stars.

###  3.1 Morphology

We start the discussion of our results with a brief comparison of the
morphological structure of the resolved ISM in the simulations with and
without runaway stars. As summarized in Table 1 the models WLM-RunP, WLM-RunM
and WLM-inplane include a treatment for runaway stars, while the simulation
WLM-fiducial is a control run that includes the feedback from massive stars
(supernovae, photoionisation and photo electric heating), as described in
Section 2.

In Fig. 1 we show face-on projections of the stellar surface density (left),
the gas surface density (center) and the thermal pressure (right) for the runs
WLM-fiducial (first row), WLM-RunP (second row), WLM-RunM (third row) and WLM-
inplane (fourth row). Qualitatively, the runs look very similar and there are
only marginal differences in the global evolution of the stellar and the ISM-
structure of the simulated systems. This trend is generally true over the full
duration of the simulation. In the simulations that include a treatment of
runaway stars, the stellar structure appears to be less clustered than in the
cases without runaway stars and the regular stellar tails that are clearly
visible in the run WLM-fiducial appear to be less prominent in the simulations
WLM-RunP, WLM-RunM and WLM-inplane. However, the effect is marginal.
In Fig. 2 we show the same results for the simulations WLM-fiducial (first
row), WLM-RunP (second row), WLM-RunM (third row) and WLM-inplane (fourth row)
for the edge-on perspective. One can easily identify the thin stellar disk.
Again we note that the overall morphology is similar in all four cases, but
one can see a slighty “puffed up” stellar disc in the runs WLM-RunP and WLM-
RunM, which include a treatment for runaway stars. This effect is not present
in the model WLM-inplane, in which velocity kicks are only applied in
directions parallel to the galactic plane. We note that the runaway stars have
a minor impact on the morphological structure of the stars and the turbulent,
supernova driven ISM, relative to the variations expected just due to the
stochastic nature of the underlying star formation and feedback prescription.

###  3.2 Star formation rate

In the top panel of Fig. 3 we show the time evolution of the star formation
rate (SFR) for all four simulations WLM-fiducial (black), WLM-RunP (blue),
WLM-RunM, and WLM-inplane (orange). Comparable to similar simulations of the
WLM system in previous work, the SFR in our fiducial run is settling around
$2\times 10^{-4}$ M⊙ yr-1, where we note that the star formation rate drops by
roughly a factor of two due to the presence of photoionising radiation,
compared to runs that only include the feedback from SNe alone (e.g. Hu et
al., 2017; Smith et al., 2021). Generally, the star formation rate in the two
runaway models WLM-RunP and RunM do not differ significantly from the star
formation rate in both the fiducial run WLM-fiducial and the run WLM-inplane,
in which we only consider velocity kicks in the x and y plane, with the
maxwellian velocity kick model. However, there are some differences that we
would like to point out in greater detail. First, there is a difference in the
strong runaway model WLM-RunP compared to the other three runs, namely the
height of the peaks in SFR which can be larger by more than a factor of two
for specific selected times. Thus the overall structure of the star formation
history in the strong runaway model appears to be more bursty. The most likely
reason for this is that the velocity kicks that we apply in this particular
case can easily drive some of the photoionising sources out of the natal
clouds. Since the molecular clouds in our WLM analogues are rather small (a
few times $10^{3}$ M⊙ up to around a maximum of $10^{4}$ M⊙) they are easily
dispersed by the effect of photoionisation. However, if a large fraction of
massive stars are quickly removed from their natal clouds, PI-radiation will
not be able to locally quench star formation as strongly as in the reference
run WLM-fiducial. This leads to a more bursty behaviour as shown in the
studies of Hu et al. (2017) and Smith et al. (2021). This is not seen to the
same degree in the simulations WLM-RunM and WLM-inplane, which adopt the more
moderate “walkaway” case with a Maxwellian velocity distribution, where the
bulk of stars (80 per cent) receive rather small velocity kicks ($\sim 10$ km
s-1) consistent with the “walkaway” picture.

In the bottom panel of Fig. 3 we show the integrated star formation rate for
all four runs, adopting the same colour scheme for the lines. While the build
up of the total stellar mass of the system is similar, the final masses (at
t=1 Gyr) are slightly different between the runs. Generally, we find that
there can be a deviation of the total mass and average star formation rate of
up to a few per cent due to model stochasticity. This is very apparent for the
runs WLM-fiducial and WLM-inplane where we find agreement on the per cent
level. However, in the runaway models WLM-RunM and WLM-RunP we find a
significant increase of the total stellar mass by 35 and 85 per cent
respectively which is a significant increase compared to model scatter of
around a few per cent. We note that we use reproducible random numbers based
on Gadget’s random number generator for the IMF sampling routine as well as
the sampling routine for the runaway stars’ velocities, our simulations are
prone to the “butterfly-effect” recently reported on in cosmological zoom
simulations (e.g. Genel et al., 2019) and molecular cloud simulations (e.g.
Keller et al., 2019). Hence our results are only strictly reproducible on the
same machine for the same number of MPI ranks and OPENMP threads. The results
of the runs WLM-fiducial and WLM-inplane are consistent with one another
within this model stochasticity, but this is not the case for the run WLM-RunP
and WLM-RunM which both show a significant boost in the integrated star
formation history.
In Fig.4 we show the supernova environmental density for all four runs,
ranging from a density of 10-4 cm-3 to 101 cm-3. The environmental density is
computed in SPH-like fashion as a kernel weighted sum over all the neighbours
of a star particle is flagged as experiencing a SN. The sum is computed before
the ejecta are distributed. For all runs we plot a total of 3000 supernova
explosions. However, due to the slight differences in the models there is
roughly a 50 Myr range in the time where the simulations reach that value,
i.e. these distributions are not plotted at the exact same time. The
distribution of the environmental density of core collapse supernovae
explosions is very similar between the different runs. There appears to be a
slight increase in the number of SN explosions in low density environments in
the model WLM-RunP compared to the other three models. Furthermore, it is
interesting to point out that the distribution of densities for the control
run WLM-inplane is very similar to the other models. From this we can directly
deduce that in the following, all major changes that we report on are driven
by “high altitude supernovae”, that explode in gas with similar density but
above or below the plane of the disk.

###  3.3 Phase structure and volume filling fraction

Next, we investigate the changes in the density-temperature phase space and
the hot volume filling fraction of the gas. In Fig. 5, we show the time
averaged density-temperature phase space diagrams over a time span of 500 Myr,
which is equivalent to 500 snapshots and around $5\times 10^{9}$ particles,
for the model WLM-fiducial(top), WLM-RunP(centre) and WLM-RunM(bottom). A more
continuously populated hot phase can be seen in both runaway models.
Additionally, we note that in both runaway cases, there is significantly more
mass in the diffuse hot phase around 105 K. Gas in HII regions is clearly
apparent as a horizontal line at 104 K.
In Fig 6 we show the same phase diagrams, also averaged over 500 Myrs, in a
volume weighted fashion (actually we weight by effective cell size which is
V1/3), which clearly reveals the volume filling nature of the hot phase of the
ISM.

In Fig. 7, we show the time evolution of the hot volume filling fraction for
all models as a function of time. To define the hot phase, we choose a cut of
$3\times 10^{4}$ K, which has two motivations. The first one is numerical, as
this is the temperature where we transit from a non-equilibrium cooling
prescription to an equilibrium cooling prescription. The second is chosen to
be comparable with previous results of Hu et al. (2017). However, we note that
this differs from the definition of the hot phase in Steinwandel et al.
(2020). We compute the hot volume filling factor as the sum over the volume
stored in the hot phase above $3\times 10^{4}$, divided by the total volume of
the simulation. The volume is computed as:

$\displaystyle V=\frac{4}{3}\pi h_{i}^{3},$ (18)

where hi is the smoothing length of particle $i$. Alternatively, one can
compute the volume as $m_{i}/\rho_{i}$. We tested both methods and note that
there is only a minor difference between the two, which typically remains
below the per cent level.

We can see from Fig 7 that the volume filling fraction of hot gas in the WLM-
RunP (strong runaway) model is significantly higher than in the other models.
The weaker runaway model WLM-RunM has a slightly higher hot gas volume filling
fraction than the WLM-fiducial and WLM-inplane models, especially at late
times. We note that the exact value for the hot phase filling fraction is very
sensitive to the exact value of the temperature cut, and can change by a
factor of five when the cut is varied within the range of of $3\times 10^{4}$K
and $3\times 10^{5}$K. Nevertheless, we find a consistently larger volume
filling factor for the hot phase independent of which temperature cut we
adopt.

Since runaway stars induce a more efficient coupling of feedback to diffuse
gas, we might expect their inclusion to effect the cold dense ISM. We do note
a substantial difference in the dense cold gas between our runs in Fig 6. In
particular, the run WLM-RunP has a noticeable dearth of volume-weighted gas
above n$=10^{3}$ cm-3. This run also has a slightly higher star formation rate
overall (as seen in Fig. 3), which suggests that the depletion time of the
cold star forming gas is shorter when runaway stars are included. In this run,
massive stars are more often ejected out of the natal cold dense gas into
warmer and diffuse regions of the ISM. Thus, the CNM and star forming gas will
be less disturbed and more likely to continue its star forming state (e.g. can
continue to be Toomre/Jeans unstable).

We now turn to a quantification of the outflow properties and loading factors.
In Fig. 8 we show the logarithm of the outflow velocity
$u\equiv\log_{10}v_{\mathrm{out}}$ as a function of the logarithm of the sound
speed $w\equiv\log_{10}c_{\mathrm{s}}$, where the colour shows the total
outflow rate in M⊙ yr-1. The gray lines indicate surfaces of constant
Bernoulli-velocity given via:

$\displaystyle
v_{\mathcal{B}}=\left(v_{\mathrm{out}}^{2}+\frac{2\gamma}{\gamma-1}c_{\mathrm{s}}^{2}\right)^{\frac{1}{2}},$
(19)

with $\gamma=5/3$. We note that these are identical to the two dimensional
PDFs presented by Kim et al. (2020) when reduced to dimensionless units in
loading factor space. Studying these diagrams for the WLM-fiducial (top) and
WLM-RunP runs reveals that the bulk of the mass is transported in the cold
phase. This is similar for both runs. However, the run WLM-RunP reveals an
extended tail in a warmer, diffuse phase. This demonstrates that runaway stars
have a clear impact on the phase structure of the resulting outflows in our
dwarf galaxy simulation.
In Fig. 9 we show the same plots for the energy outflow rate, for WLM-fiducial
on the top and WLM-RunP on the bottom, indicating that while there is a cold
energy transporting phase in our dwarf galaxy simulation, there is a tail that
extends to higher specific energy. This is even more true in the bottom panel
for strong runaway case WLM-RunP, where we find a tail that not only extends
to higher sound speed but also to higher outflow velocity, which is direct
evidence that runaway stars can not only boost the mass flux but can
additionally increase the specific energy loading as well as the total energy
loading of the wind significantly.

###  3.4 Outflow rate and outflow mass loading

We plot the main outflow diagnostics for all four simulations in Fig. 10. The
top left panel shows the mass outflow rate (M⊙ yr-1), the top right panel
shows the metal outflow rate (M⊙ yr-1), the bottom left panel shows the
momentum outflow rate (M⊙ km s-1 yr-1) and the bottom left panel shows the
energy outflow rate (erg yr-1). We note that the outflow rates are measured at
a height of $z=1$ kpc above the midplane in a patch of thickness 100 pc. We
find a slight boost in the mass and the metal outflow rate by a factor of
around two in the late evolutionary stages of the galactic disc when runaway
stars are included. The effect is roughly the same between the strong runaway
(WLM-RunP) and weaker runaway (WLM-RunM) case. The trend is slightly stronger
in the metal outflow rate, which suggests that that the runaway stars deposit
their ejecta directly into the outflowing gas. While the mass and metal
outflow rates are already clearly boosted in the runs with runaway stars, the
effect is even stronger for the momentum and energy outflow rates. For the
momentum outflow rate we find a boost of around a factor of five, while we see
a boost in the energy outflow rate by at least a factor of five and in single
peaks up to a factor of around ten. Comparing with the reference runs WLM-
fiducial and WLM-inplane clearly indicates that this boost stems from “high-
altitude supernovae” in low density gas around the disc scale height of the
gas, as this is the mean travel distance of the fast moving runaway stars. The
effect is slightly suppressed in the WLM-RunM model which is consistent with
the scenario of walkaway stars.

In Fig.11 we show the mass (top left), metal (top right), momentum (bottom
left) and energy (bottom right) loading for all four simulations considered in
this paper. While we see similar trends for the loading factors, we note that
the trend in the mass loading is slightly suppressed compared to the boost we
find in the mass outflow rate.

###  3.5 Inflow rates

In Fig.12 we show the inflow rate for mass (top left), metals (top right),
momentum (bottom left) and energy (bottom right) for the simulations WLM-
fiducial (black), WLM-RunP (blue), WLM-RunM, and WLM-inplane (orange). We
measure the inflow rate in the same 100 pc thick patch located at 1 kpc
above/below the midplane. Generally, one can see that in all the models, the
inflow rates decrease with time, then increase slightly after about 500 Myr.
This increase is more pronounced in the models with runaway stars. In the
models with runaway stars, the mass, metal and momentum inflow rates show a
similar degree of increase as the increase in outflow rates. However, the
energy inflow rate decreases by only a factor of 2-3, not a factor of $\sim
10$ as was seen in the energy outflow rate for the strong runaway model. This
implies that the SN-feedback energy of “high altitude supernovae” is either
efficiently kept in the CGM due to mixing with lower specific energy gas, or
lost to cooling in the highly metal rich outflowing material. The latter case
seems unlikely when we consider the phase-structure of the mass and volume
weighted density-temperature diagrams. These reveal more mass in the diffuse
105 K regime at low densities, indicating that there is an slight increase in
the virial temperature of the system when runaway stars are included. We note
that inflow rates are slightly larger in the runs that include a
runaway/walkaway model. This can be understood in combination with the
slightly increased mass outflow rates we observe in the runaway models. We
find that the peaks in outflow and inflow rate are shifted by roughly 150-200
Myr, indicating that some fraction of the gas is actually “recylced” over a
few dynamical times of our simulated dwarf.

###  3.6 Metal enrichment factor

Finally, we discuss the metal outflow enrichment factor $y_{Z}^{\mathrm{out}}$
and the metal inflow enrichment factor $y_{Z}^{\mathrm{in}}$, which we show
for all simulations in Fig.13 on the top and the bottom respectively. There is
a striking increase of the metal outflow enrichment factor when runaway star
modelling is included in the simulations. Generally we find a boost of around
15 to 20 percent in the metal outflow enrichment factor compared to the
fiducial run and the run where we only apply the velocity kicks in the x- and
y-plane. We note that while the metal inflow enrichment factor is also boosted
compared to the reference simulations the effect is not as strong as for the
outflow enrichment factor. This indicates that runaway stars can significantly
contribute to the metal enrichment of the CGM and can boost its metallicity by
about 15 percent. Since metals are the most important coolant in the low
density regime, this might have crucial consequences for mixing properties in
the CGM. However, since we do not explicitly model the CGM in our simulation
by including a hot cooling halo, we postpone a detailed study of CGM cooling,
heating and mixing rates alongside with a detailed analysis of metal emission
line properties to future work.

Figure 13: We show the outflow metal enrichment factor $y_{Z}^{\mathrm{out}}$
(top) and the inflow metal enrichment factor $y_{Z}^{\mathrm{in}}$ for the
runs WLM-fiducial (black), WLM-RunP (blue), WLM-RunM (red) and WLM-inplane
(orange). We find a significant boost in the outflow metal enrichment factor
in the runs that include runaway stars. Naturally, the inflowing gas is more
metal enriched as well compared to the background ISM metallicity.

##  4 Discussion

###  4.1 Interpretation and comparison to previous work

We start our discsussion by briefly comparing to previous results of simialr
WLM systems without the inclusion of runaway star modeling, which has been
used in several studies throughout the literature (Hu et al., 2016, 2017;
Emerick et al., 2019; Hu et al., 2019; Gutcke et al., 2021; Smith et al.,
2021; Hislop et al., 2022). The most appropriate of these to compare to are
the simulations by Hu et al. (2017) and Smith et al. (2021) since the
implemented physics modules cover the same range from the cooling and SN-
feedback routines to the prescription for photo-electric heating and
photoionisation. This is important since Hu et al. (2017) used a pressure-
energy SPH formalism with 100 neighbours and Smith et al. (2021) adopted the
moving mesh code arepo for their simulations which employs a second order
Riemann-method on an unstaggered moving mesh. The latter is numerically more
similar to our MFM approach. Generally, all codes agree quite well for our
fiducial model when compared to the appropriate runs of Hu et al. (2017) (PI-
SN-PE) and Smith et al. (2021) (low$\Sigma$ in their Appendix A). This is true
across star formation histories and loading factors. A more detailed study on
the subtle differences between SPH and MFM on identical setups will be
investigated in future work.

Our results indicate that even conservative modelling of runaway stars in
galaxy simulations with a resolved treatment of the ISM and stellar feedback
can have a significant impact on the properties of large scale stellar driven
winds. This most apparent in the energy outflow and energy loading factors,
which we show in our respective Fig. 10 and Fig. 11. It is remarkable to point
out that the inclusion of runaway stars can boost the energy loading factor by
up to an order of magnitude, up to a value of $\sim 0.1$ in the peak values.
On average, we still find an increase by roughly a factor of 5 in the model
WLM-runP and a factor of 3 in the model WLM-runM. We note that previous work,
such as from the tigress suite (e.g. Kim & Ostriker, 2015, 2018; Kim et al.,
2020) or from (e.g. Li et al., 2017), typically find a value of the energy
loading of around $\sim 0.1$ for Milky Way like conditions. Our results
suggest that similarly high values of the energy loading might also be found
in dwarf galaxies. The imprints of the runaway star process are not very clear
in the mass and temporal weighted density-temperature phase space diagrams
that we show in Fig.5, but become more apparent in the volume weighted phase-
space diagrams shown in Fig.6. There we can clearly see that the volume of the
hot phase significantly increases when any runaway modelling is included. This
can be quantified not only as an increase in the specific energy loading but
also in the total energy loading, as we can see in Fig.9. While the effect on
the mass loading is somewhat weaker, we still see a significant increase in
the mass flux rates when runaway stars are included (Fig. 8).

In this context it is interesting to point out that while we do not explicitly
include a CGM around our dwarf galaxy simulation, it is not that case that
there is no pristine low density gas around the dwarf due to the 500 Myr of
initial evolution with turbulent SN-driving that we apply before examining any
science results. This can be clearly deduced from Fig. 4 where we find little
difference in the SN-environmental distribution between the fiducial run WLM-
fiducial and the run WLM-inplane where the runaway kicks are only applied in
the plane of the galaxy. The similar shape of these implies two things
directly. First, the increase in the specific energy of the wind is not driven
by the re-positioning of the SNe in the disc due to the presence of the
runaway stars but rather stems from the “high-altitude” SNe, introduced by the
presence of the runaway star modelling. That this effect is still present in
the weaker walkaway case is a rather promising result as it indicates that the
re-positioning of SNe in the vertical direction is more important at constant
environmental SN-density. Second, the fact that the low density distribution
of SN explosion environment is only weakly affected by the presence of any
runaway modelling indicates that the SNe have to explode around the disc scale
height in gas that is still mostly resolved even without an additional CGM
that would provide an external pressure. If this were not the case, the
distribution of environmental SN-densities would be significantly altered and
shift to an overabundance in events in unresolved low density gas.
Finally, it is important to point out that we find a significant increase in
the metal mass loading, by around 20 per cent in both runaway models as we
show in 13. This can have crucial consequences for the cooling behaviour
around dwarf galaxy systems and observational properties of the CGM.
We can compare our results with two notable studies concerning runaway stars.
One was presented by Andersson et al. (2020), and was carried out with the
adaptive mesh refinement code Ramses (Teyssier, 2002) for a global MW-
simulation. Another was carried out by Kim & Ostriker (2018) within the
tigress framework of local stratified environments. At face value, these
studies seem to disagree with one another, as Kim & Ostriker (2018) reports
that runaways have a very weak effect, while Andersson et al. (2020) reports a
significant boost in the wind outflow rates when runaway modeling is included.
However, the comparison is tricky due to the different setup choices.
Furthermore, both studies are evolved for a short amount of time compared to
any global MW-scale. The comparison to our work is even more complicated since
our simulations are designed to study very low mass dwarf galaxies at very
high resolution. However, Andersson et al. (2020) state that they find that
runaway stars do not have any effect in dwarf galaxies, which they discuss
very briefly without showing the details of the dwarf simulations. These
appear to be centred around very low mass cosmological dwarf galaxy
simulations taken from the edge sample (Agertz et al., 2020). Indeed, we also
find that if we apply the runaway modelling in lower mass isolated discs we do
not see a strong effect on any of the properties, which may be because these
do not form a thin disc due to the strong impact of SN-feedback. Hence our
results are not necessarily in disagreement with what Andersson et al. (2020)
report, but we represent difference in the galaxy properties between our study
and theirs. Furthermore, we note that some of the dwarf galaxy simulations
that are carried out with ramses by now do not see an effect of runway stars
when only SN-feedback is employed (Andersson p.c. and Andersson et al. in
prep.) The same can be said about a comparison with Kim & Ostriker (2018),
where again, the results are very hard to compare since we are in very
different outflow regime, as can be seen from the shape of our Fig. 8 and Fig.
9 compared to their results as shown in Kim et al. (2020). For now we can only
conclude that our results seem to be somewhere in between the studies of
Andersson et al. (2020) and Kim & Ostriker (2018), with runaways having a
weaker impact than in the Milky Way simulations of Andersson et al. (2020),
but stronger than in the work of Kim & Ostriker (2018). However, differently
from both studies we find a significant increase in the appearance of high
specific energy outflows although we note that while Kim & Ostriker (2018) is
investigating energy loading as well, this is not put forward in Andersson et
al. (2020). Nevertheless, we recommend that future resolved simulations
include some treatment of runaways to capture this effect. However, we note
that even with the runaway star modelling we do not achieve the high loading
factors observed in the supernova only case (e.g. Hu et al., 2016, 2017). The
lower loading factors in out case are mostly driven by the inclusion of
photoionising radiation and the runway star implementation seem to be able to
counterbalance that effect on the loading factors.

###  4.2 Model limitations

In this sub-section, we briefly discuss some caveats and model limitations of
our simulations. First, our two presented runaway models strictly only
represent a natal kicks in a binary star formation scenario. However, we note
that we only sample a single star IMF that strictly speaking does not account
for a fraction of the stars to be formed in binary or few body systems. The
latter would allow us to study more complex runaway scenarios like ejection of
massive stars due to few body interactions or ejection of B stars in a OB
binary due to supernova kicks when the O-star explodes. However, this could be
realised by adding a time delay to the kick, which is straightforward to
implement in our code. It is not very likely that this will have a major
impact on the conclusions of our study.
Additionally, we do not include an explicit pre-existing CGM around our dwarf
galaxy, which could result in an overestimate of the boost in the outflow
rates and loading factors. However, there are a few arguments that support our
picture even without the explicit inclusion of a low density CGM around the
dwarf. First, dwarfs have a very low virial temperature of $\sim 10^{4}$ K. In
combination with the assumption of a very low density in the CGM compared to
the galactic disc, it is rather unlikely that the CGM of dwarf galaxies is
heavily supported by thermal pressure which would work against the thermal
pressure provided by high altitude SNe due to the presence of runaway stars.
Second, we note that the SN-environmental density distributions in all runs
are very similar and there is no excess of SNe in lower density regions. This
directly indicates that even above the disc scale height the resolution seems
still to be sufficient to model the evolution of individual SN-remnants.

##  5 Conclusions and Outlook

###  5.1 Summary

Table 2: Summary of the mean outflow diagnostics once the system has reached a quasi-self-regulating state Name | SFR [M⊙ yr-1] |  $\dot{M}_{\mathrm{out}}$ [M⊙ yr-1] |  $\dot{M}_{\mathrm{Z,out}}$ [M⊙ yr-1] |  $\dot{p}_{\mathrm{out}}$ [M⊙ km s-1 yr-1] |  $\dot{E}_{\mathrm{out}}$ [$10^{43}$ erg yr-1]
---|---|---|---|---|---
WLM-fiducial | 0.00023 | 0.00632 | 0.00068 | 0.194 | 2.225
WLM-RunP | 0.00021 | 0.00816 | 0.00100 | 0.291 | 15.298
WLM-RunM | 0.00029 | 0.00757 | 0.00087 | 0.243 | 8.512
WLM-inplane | 0.00024 | 0.00686 | 0.00074 | 0.200 | 3.085
Table 3: Summary of the mean loading factors once the system has reached a quasi-self-regulating state Name | $\eta_{\mathrm{M}}^{\mathrm{out}}$ | $\eta_{\mathrm{Z}}^{\mathrm{out}}$ | $\eta_{\mathrm{p}}^{\mathrm{out}}$ | $\eta_{\mathrm{E}}^{\mathrm{out}}$ | $y_{\mathrm{Z}}^{\mathrm{out}}$
---|---|---|---|---|---
WLM-fiducial | 31.4 | 125.2 | 2.93 | 0.11 | 1.089
WLM-RunP | 45.7 | 136.4 | 3.28 | 0.5 | 1.23
WLM-RunM | 39.0 | 143.5 | 3.34 | 0.35 | 1.15
WLM-inplane | 34.0 | 133.0 | 2.97 | 0.13 | 1.088

We carried out high resolution numerical simulations of isolated dwarf
galaxies with and without a treatment for runaway stars to probe the effect of
this process on galactic outflows and related galaxy and ISM properties. We
showed that the effect of runaway stars is weak with respect to the resulting
structure of the multiphase ISM. Furthermore, we showed that there is an
impact on the global star formation rate and the build up of the total stellar
mass of the system as a function of time for all simulations that include
runaway stars compared to the fiducial run. While this results in very similar
distributions for the PDFs of the environmental density where SNe explode, the
run WLM-RunP shows a slight excess of supernovae in lower density
environments, which is beyond the intrinsic model scatter. However, the fact
that all the models, including the run WLM-inplane, show similar SN density
distributions is clear evidence that the main changes in global outflow
evolution due to runaways stars is caused by the stars that explode at higher
altitude around the midplane. The fact that this trend remains for the weaker
runaway model WLM-RunM indicates that runaway stars can heat the gas at the
important boundary between galactic disc and CGM.
This can be quantified by investigating the hot phase of the ISM in density-
temperature phase space in the cases with and without runaway stars. Averaging
over about half of the simulation run time, we find an excess in the mass of
hot gas in the diffuse low density ISM in the simulations with runaway stars.
In Fig. 10 we show the mass, metal, momentum and energy outflow rates as a
function of time. All quantities are measured at a height of 1 kpc above the
midplane in a thin slice with thickness 0.1 kpc. We find a slight boost in the
overall mass outflow rate which is slightly more dominant when only the metal
mass is considered, showing first evidence that runaway stars can lead to an
excess of metals in the outflowing gas. The more dominant effect of the
runaway stars can be seen in the momentum and energy outflow rates, where the
momentum outflow rate is boosted by at least a factor of two in the models
WLM-RunP and WLM-RunM. The effect is even more dramatic for the energy outflow
rate where in the model WLM-RunP it is boosted by almost one order of
magnitude and in the model WLM-RunM it increases by roughly a factor of 7. As
this increase is not seen in the fiducial run WLM-fiducial and the run WLM-
inplane in which we apply the natal kicks only in the midplane of the galaxy,
this provides clear evidence that runaway stars could contribute to the
driving of high specific energy winds. We note that while these trends are
also seen in the momentum and energy loading they are somewhat reduced for the
mass and metal outflow loading. There is a clear trend of stronger metal
enrichment in the models that include a treatment of runaway stars. We
summarise the mean outflow parameters after saturation in Table 2 and the mean
outflow loading factors in Table 3. The key findings of this paper can be
summarised as follows

  1. 1.

Runaway stars affect the structure of the multiphase turbulent ISM and can
lead to a slightly “puffier” stellar disc.

  2. 2.

The global star formation rate and the distribution of SN-environmental
densities is increasing in the simulations with runaway stars.

  3. 3.

The cold ISM gas is less disturbed when runaways are included. This leads to
enhanced star formation as the cold gas can continue to be gravitationally
unstable.

  4. 4.

Runaways stars can contribute to the build-up of the hot phase of the ISM and
promote outflow launching.

  5. 5.

Runaway stars can boost the mass outflow rate and metal outflow rate while
keeping mass outflow and metal outflow loading roughly constant compared to
the fiducial run WLM-fiducial

  6. 6.

Runaways stars can boost the momentum outflow rate by around a factor of at
least 2, a trend that remains for the momentum outflow loading.

  7. 7.

Runaways stars can boost the energy outflow rate by around a factor of 7-10 in
the peaks and around 3-5 in the mean, a trend that remains for the momentum
outflow loading. Therefore, they can potentially contribute to the origin of
high specific energy winds in dwarf galaxies.

###  5.2 Outlook and future work

We close by briefly discussing a few interesting ideas and implications for
future work. While we find a rather weak impact of runaway stars on the mass
outflow rate and specifically the outflow mass loading, runaway stars could
have implications beyond what is represented in our study. Recently, some
groups have explored cosmic ray (CRs) driven winds (e.g. Pakmor et al., 2016;
Buck et al., 2020; Hopkins et al., 2020, 2021c, 2021d, 2021e, 2021a; Hopkins
et al., 2021b; Hopkins et al., 2022; Hopkins, 2022). In particular, Quataert
et al. (2022b) and Quataert et al. (2022a) study in detail CR wind driving
scenarios for CR-streaming and CR-diffusion. The outcome of this extensive
analytic study is that winds by CR-streaming require a significant Alfvén
velocity $v_{\mathrm{A}}$ and thus magnetic field strength. This is
intrinsically hard to obtain in dwarf galaxies as the ISM turbulence is rather
weak due to low SN-rates, which limits small-scale turbulent dynamo action.
Furthermore, since dwarfs are often only weakly rotationally supported, there
is no large scale dynamo possible. On the other hand, the diffusion driven
winds require a large diffusion coefficient in order to obtain consistency
with gamma-ray observations from pion decays that originate from cosmic ray
protons, which puts rather tight constraints on the exact value of this
coefficient. An interesting direction to explore in future work (both
numerically and analytically) would be to constrain the effect of cosmic rays
seeded by runaway stars at or above the wind launching scale which could
significantly contribute to outflows driven by cosmic rays above the midplane.
Hence, one could imagine runaways stars as an external boost for the cosmic
ray diffusion coefficient. This could be explored with future simulations and
improved analytic modelling.

## Data Availability Statement

The data used in this article will be made available based on reasonable
request to the corresponding author.

## Acknowledgements

UPS is supported by a Flatiron Research Fellowship (FRF) at the Center of
Computational Astrophysics at the Flatiron Institute. The Flatiron Institute
is supported by the Simons Foundation. RSS acknowledges support from the
Simons Foundation. B.B. is grateful for funding support from the Simons
Foundation, Sloan Foundation, and the Packard Foundation. UPS acknowledges the
useful discussions with Eric Andersson, Drummond Fielding, Chang-Goo Kim,
Thorsten Naab, Eve Ostriker and Romain Teyssier. UPS acknowledges the
computing time provided by the Leibniz-Rechenzentrum (LRZ) in Garching on
SuperMUC-NG under the project number pn72bu, as well as the computing time on
the c2pap-cluster in Garching under the project number pr27mi on which most
the simulations presented here have been carried out.

## References

  * Agertz et al. (2013) Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, ApJ, 770, 25
  * Agertz et al. (2020) Agertz O., et al., 2020, MNRAS, 491, 1656
  * Andersson et al. (2020) Andersson E. P., Agertz O., Renaud F., 2020, MNRAS, 494, 3328
  * Bakes & Tielens (1994) Bakes E. L. O., Tielens A. G. G. M., 1994, ApJ, 427, 822
  * Bergin et al. (2004) Bergin E. A., Hartmann L. W., Raymond J. C., Ballesteros-Paredes J., 2004, ApJ, 612, 921
  * Blaauw (1961) Blaauw A., 1961, Bull. Astron. Inst. Netherlands, 15, 265
  * Buck et al. (2020) Buck T., Pfrommer C., Pakmor R., Grand R. J. J., Springel V., 2020, MNRAS, 497, 1712
  * Chieffi & Limongi (2004) Chieffi A., Limongi M., 2004, ApJ, 608, 405
  * Comerón & Pasquali (2007) Comerón F., Pasquali A., 2007, A&A, 467, L23
  * Eldridge et al. (2011) Eldridge J. J., Langer N., Tout C. A., 2011, MNRAS, 414, 3501
  * Emerick et al. (2019) Emerick A., Bryan G. L., Mac Low M.-M., 2019, MNRAS, 482, 1304
  * Fielding et al. (2017) Fielding D., Quataert E., Martizzi D., Faucher-Giguère C.-A., 2017, MNRAS, 470, L39
  * Forbes et al. (2016) Forbes J. C., Krumholz M. R., Goldbaum N. J., Dekel A., 2016, Nature, 535, 523
  * Fujii & Portegies Zwart (2011) Fujii M. S., Portegies Zwart S., 2011, Science, 334, 1380
  * Fujii et al. (2021) Fujii M. S., Saitoh T. R., Hirai Y., Wang L., 2021, PASJ, 73, 1074
  * Gaburov & Nitadori (2011) Gaburov E., Nitadori K., 2011, MNRAS, 414, 129
  * Gatto et al. (2017) Gatto A., et al., 2017, MNRAS, 466, 1903
  * Genel et al. (2019) Genel S., et al., 2019, ApJ, 871, 21
  * Georgy et al. (2013) Georgy C., et al., 2013, A&A, 558, A103
  * Gies (1987) Gies D. R., 1987, ApJS, 64, 545
  * Girichidis et al. (2016) Girichidis P., et al., 2016, MNRAS, 456, 3432
  * Girichidis et al. (2022) Girichidis P., Pfrommer C., Pakmor R., Springel V., 2022, MNRAS, 510, 3917
  * Glover & Clark (2012) Glover S. C. O., Clark P. C., 2012, MNRAS, 421, 9
  * Glover & Mac Low (2007) Glover S. C. O., Mac Low M.-M., 2007, ApJS, 169, 239
  * Grond et al. (2019) Grond J. J., Woods R. M., Wadsley J. W., Couchman H. M. P., 2019, MNRAS, 485, 3681
  * Grudić et al. (2022) Grudić M. Y., Guszejnov D., Offner S. S. R., Rosen A. L., Raju A. N., Faucher-Giguère C.-A., Hopkins P. F., 2022, arXiv e-prints, p. arXiv:2201.00882
  * Guedes et al. (2011) Guedes J., Callegari S., Madau P., Mayer L., 2011, ApJ, 742, 76
  * Guszejnov et al. (2022) Guszejnov D., Markey C., Offner S. S. R., Grudić M. Y., Faucher-Giguère C.-A., Rosen A. L., Hopkins P. F., 2022, arXiv e-prints, p. arXiv:2201.01781
  * Gutcke et al. (2021) Gutcke T. A., Pakmor R., Naab T., Springel V., 2021, MNRAS, 501, 5597
  * Haid et al. (2016) Haid S., Walch S., Naab T., Seifried D., Mackey J., Gatto A., 2016, MNRAS, 460, 2962
  * Haid et al. (2018) Haid S., Walch S., Seifried D., Wünsch R., Dinnbier F., Naab T., 2018, MNRAS, 478, 4799
  * Hirai et al. (2021) Hirai Y., Fujii M. S., Saitoh T. R., 2021, PASJ, 73, 1036
  * Hirschmann et al. (2014) Hirschmann M., Dolag K., Saro A., Bachmann L., Borgani S., Burkert A., 2014, MNRAS, 442, 2304
  * Hislop et al. (2022) Hislop J. M., Naab T., Steinwandel U. P., Lahén N., Irodotou D., Johansson P. H., Walch S., 2022, MNRAS, 509, 5938
  * Hoogerwerf et al. (2000) Hoogerwerf R., de Bruijne J. H. J., de Zeeuw P. T., 2000, ApJ, 544, L133
  * Hoogerwerf et al. (2001) Hoogerwerf R., de Bruijne J. H. J., de Zeeuw P. T., 2001, A&A, 365, 49
  * Hopkins (2015) Hopkins P. F., 2015, MNRAS, 450, 53
  * Hopkins (2022) Hopkins P. F., 2022, arXiv e-prints, p. arXiv:2202.05283
  * Hopkins et al. (2012) Hopkins P. F., Quataert E., Murray N., 2012, MNRAS, 421, 3522
  * Hopkins et al. (2014) Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 445, 581
  * Hopkins et al. (2018) Hopkins P. F., et al., 2018, MNRAS, 480, 800
  * Hopkins et al. (2020) Hopkins P. F., et al., 2020, MNRAS, 492, 3465
  * Hopkins et al. (2021a) Hopkins P. F., Butsky I. S., Panopoulou G. V., Ji S., Quataert E., Faucher-Giguere C.-A., Keres D., 2021a, arXiv e-prints, p. arXiv:2109.09762
  * Hopkins et al. (2021b) Hopkins P. F., Squire J., Butsky I. S., Ji S., 2021b, arXiv e-prints, p. arXiv:2112.02153
  * Hopkins et al. (2021c) Hopkins P. F., Chan T. K., Ji S., Hummels C. B., Kereš D., Quataert E., Faucher-Giguère C.-A., 2021c, MNRAS, 501, 3640
  * Hopkins et al. (2021d) Hopkins P. F., Chan T. K., Squire J., Quataert E., Ji S., Kereš D., Faucher-Giguère C.-A., 2021d, MNRAS, 501, 3663
  * Hopkins et al. (2021e) Hopkins P. F., Squire J., Chan T. K., Quataert E., Ji S., Kereš D., Faucher-Giguère C.-A., 2021e, MNRAS, 501, 4184
  * Hopkins et al. (2022) Hopkins P. F., Squire J., Butsky I. S., 2022, MNRAS, 509, 3779
  * Hu et al. (2014) Hu C.-Y., Naab T., Walch S., Moster B. P., Oser L., 2014, preprint,  (arXiv:1402.1788)
  * Hu et al. (2016) Hu C.-Y., Naab T., Walch S., Glover S. C. O., Clark P. C., 2016, MNRAS, 458, 3528
  * Hu et al. (2017) Hu C.-Y., Naab T., Glover S. C. O., Walch S., Clark P. C., 2017, MNRAS, 471, 2151
  * Hu et al. (2019) Hu C.-Y., Zhukovska S., Somerville R. S., Naab T., 2019, arXiv e-prints,
  * Hu et al. (2021) Hu C.-Y., Sternberg A., van Dishoeck E. F., 2021, ApJ, 920, 44
  * Hu et al. (2022) Hu C.-Y., Schruba A., Sternberg A., van Dishoeck E. F., 2022, arXiv e-prints, p. arXiv:2201.03885
  * Keller et al. (2019) Keller B. W., Wadsley J. W., Wang L., Kruijssen J. M. D., 2019, MNRAS, 482, 2244
  * Kim & Ostriker (2015) Kim C.-G., Ostriker E. C., 2015, ApJ, 802, 99
  * Kim & Ostriker (2018) Kim C.-G., Ostriker E. C., 2018, ApJ, 853, 173
  * Kim et al. (2017) Kim C.-G., Ostriker E. C., Raileanu R., 2017, ApJ, 834, 25
  * Kim et al. (2019) Kim J.-G., Kim W.-T., Ostriker E. C., 2019, ApJ, 883, 102
  * Kim et al. (2020) Kim C.-G., et al., 2020, ApJ, 903, L34
  * Kroupa (1995) Kroupa P., 1995, MNRAS, 277, 1491
  * Lahén et al. (2019a) Lahén N., Naab T., Johansson P. H., Elmegreen B., Hu C.-Y., Walch S., Steinwand el U. P., Moster B. P., 2019a, arXiv e-prints, p. arXiv:1911.05093
  * Lahén et al. (2019b) Lahén N., Naab T., Johansson P. H., Elmegreen B., Hu C.-Y., Walch S., 2019b, ApJ, 879, L18
  * Lancaster et al. (2021a) Lancaster L., Ostriker E. C., Kim J.-G., Kim C.-G., 2021a, ApJ, 914, 89
  * Lancaster et al. (2021b) Lancaster L., Ostriker E. C., Kim J.-G., Kim C.-G., 2021b, ApJ, 914, 89
  * Lancaster et al. (2021c) Lancaster L., Ostriker E. C., Kim J.-G., Kim C.-G., 2021c, ApJ, 922, L3
  * Lejeune et al. (1997) Lejeune T., Cuisinier F., Buser R., 1997, A&AS, 125, 229
  * Lejeune et al. (1998) Lejeune T., Cuisinier F., Buser R., 1998, A&AS, 130, 65
  * Li et al. (2017) Li M., Bryan G. L., Ostriker J. P., 2017, ApJ, 841, 101
  * Malkov & Zinnecker (2001) Malkov O., Zinnecker H., 2001, MNRAS, 321, 149
  * Marinacci et al. (2014) Marinacci F., Pakmor R., Springel V., 2014, MNRAS, 437, 1750
  * Martizzi et al. (2015) Martizzi D., Faucher-Giguère C.-A., Quataert E., 2015, MNRAS, 450, 504
  * Naab & Ostriker (2017) Naab T., Ostriker J. P., 2017, ARA&A, 55, 59
  * Nelson et al. (2018) Nelson D., et al., 2018, MNRAS, 475, 624
  * Nelson et al. (2019) Nelson D., et al., 2019, MNRAS, 490, 3234
  * Oh & Kroupa (2016) Oh S., Kroupa P., 2016, A&A, 590, A107
  * Orr et al. (2021) Orr M. E., Fielding D. B., Hayward C. C., Burkhart B., 2021, arXiv e-prints, p. arXiv:2109.14656
  * Orr et al. (2022) Orr M. E., Fielding D. B., Hayward C. C., Burkhart B., 2022, ApJ, 924, L28
  * Pakmor et al. (2016) Pakmor R., Pfrommer C., Simpson C. M., Springel V., 2016, ApJ, 824, L30
  * Perets & Šubr (2012) Perets H. B., Šubr L., 2012, ApJ, 751, 133
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 473, 4077
  * Quataert et al. (2022a) Quataert E., Jiang F., Thompson T. A., 2022a, MNRAS, 510, 920
  * Quataert et al. (2022b) Quataert E., Thompson T. A., Jiang Y.-F., 2022b, MNRAS, 510, 1184
  * Rantala et al. (2017) Rantala A., Pihajoki P., Johansson P. H., Naab T., Lahén N., Sawala T., 2017, ApJ, 840, 53
  * Rantala et al. (2020) Rantala A., Pihajoki P., Mannerkoski M., Johansson P. H., Naab T., 2020, MNRAS, 492, 4131
  * Renzo et al. (2019) Renzo M., et al., 2019, A&A, 624, A66
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Silva & Napiwotzki (2011) Silva M. D. V., Napiwotzki R., 2011, MNRAS, 411, 2596
  * Smith et al. (2021) Smith M. C., Bryan G. L., Somerville R. S., Hu C.-Y., Teyssier R., Burkhart B., Hernquist L., 2021, MNRAS, 506, 3882
  * Somerville & Davé (2015) Somerville R. S., Davé R., 2015, ARA&A, 53, 51
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Springel et al. (2005) Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776
  * Steinwandel et al. (2020) Steinwandel U. P., Moster B. P., Naab T., Hu C.-Y., Walch S., 2020, MNRAS,
  * Teyssier (2002) Teyssier R., 2002, A&A, 385, 337
  * Vogelsberger et al. (2014) Vogelsberger M., et al., 2014, MNRAS, 444, 1518
  * Walch et al. (2015) Walch S., et al., 2015, MNRAS, 454, 238
  * Westera et al. (2002) Westera P., Lejeune T., Buser R., Cuisinier F., Bruzual G., 2002, A&A, 381, 524
  * Wolfire et al. (2003) Wolfire M. G., McKee C. F., Hollenbach D., Tielens A. G. G. M., 2003, ApJ, 587, 278