# How non-thermal pressure impacts the modelling of star formation in galaxy
formation simulations

Eirini Batziou,1,2 Ulrich P. Steinwandel,3 Klaus Dolag2,1 and Milena
Valentini4,5,6,2
1Max Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1, 85748,
Garching, Germany
2Universitäts-Sternwarte München, Fakultät für Physik, LMU Munich,
Scheinerstr. 1, 81679, Germany
3Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue,
New York, NY 10010, USA
4Astronomy Unit, Department of Physics, University of Trieste, via Tiepolo 11,
I-34131 Trieste, Italy
5INAF - Osservatorio Astronomico di Trieste, via Tiepolo 11, I-34131 Trieste,
Italy
6ICSC - Italian Research Center on High Performance Computing, Big Data and
Quantum Computing E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

In cosmological simulations of large-scale structure star formation and
feedback in galaxies are modelled by so-called sub-grid models, that represent
a physically motivated approximation of processes occurring below the
resolution limit. However, when additional physical processes are considered
in these simulations, for instance, magnetic fields or cosmic rays, they are
often not consistently coupled within the descriptions of the underlying sub-
grid star formation models. Here, we present a careful study on how one of the
most commonly used sub-grid models for star formation in current large-scale
cosmological simulations can be modified to self consistently include the
effects of non-thermal components (e.g., magnetic fields) within the fluid. We
demonstrate that our new modelling approach, that includes the magnetic
pressure as an additional regulation on star formation, can reproduce global
properties of the magnetic field within galaxies in a setup of an isolated
Milky Way-like galaxy simulation, but is also successful in reproducing local
properties such as the anti-correlation between the local magnetic field
strength with the local star formation rate as observed in galaxies (i.e. NGC
1097). This reveals how crucial a consistent treatment of different physical
processes is within cosmological simulations and gives guidance for future
simulations.

###### keywords:

galaxies – galaxies: magnetic fields – galaxies: star formation – numerical
simulations

††pubyear: 2022††pagerange: How non-thermal pressure impacts the modelling of
star formation in galaxy formation simulations–How non-thermal pressure
impacts the modelling of star formation in galaxy formation simulations

##  1 Introduction

Numerical simulations of structure formation are essential for understanding
the dynamics of physical processes in the Universe. Usually, they describe
scales of Giga-parsec (Gpc) and include galaxy clusters, the intracluster
medium (ICM), groups of galaxies, galaxies and their constituents. Modelling
different physical processes in such a variety of spatial and mass scales is a
challenge that brings the necessity to use physically motivated models
describing the physics on scales that are not possible to resolve given the
current computational resources.

Understanding galaxy formation and evolution has been a long-standing problem,
theoretically approached by semi-analytic methods (e.g. Lacey & Silk, 1991;
Kauffmann et al., 1994), hydrodynamical simulations of galaxy groups (e.g.,
Hirschmann et al., 2014; Vogelsberger et al., 2014; Schaye et al., 2015;
Pillepich et al., 2018; Crain & van de Voort, 2023), and hydrodynamical
simulations of single galaxies (e.g., Governato et al., 2008; Wang et al.,
2015; Grand et al., 2017; Hopkins et al., 2014; Oser et al., 2010; Guedes et
al., 2011; Aumer et al., 2013; Marinacci et al., 2014; Murante et al., 2015;
Agertz & Kravtsov, 2016; Valentini et al., 2017, 2018; Valentini et al., 2023;
Hopkins, 2015; Giammaria et al., 2021; Hopkins et al., 2018, 2022, 2023).

Magnetic fields are ubiquitous in the Universe, and observations of disk
galaxies show the presence of magnetic fields in (rough) energy equipartition
with the remaining energy components of the interstellar medium (ISM) (see,
e.g., Beck, 2015, for a review). Therefore, there is a growing effort to
involve magnetic fields in simulations of galaxies (e.g., Wang & Abel, 2009;
Dubois & Teyssier, 2010; Pakmor & Springel, 2013; Rieder & Teyssier, 2016;
Butsky et al., 2017; Steinwandel et al., 2019, 2020; Wibking & Krumholz,
2023a) and in cosmological zoom-in simulations (e.g., Beck et al., 2012;
Rieder & Teyssier, 2017a, b; Pakmor et al., 2017; Su et al., 2017; Martin-
Alvarez et al., 2018; Hopkins et al., 2020; Pakmor et al., 2024; Martin-
Alvarez et al., 2021). Meanwhile, single large-scale cosmological simulations
have started to investigate the magnetic field structure on larger scales
(e.g., Pillepich et al., 2018; Marinacci et al., 2018; Ramesh et al., 2023).

Star formation and stellar feedback are the key ingredients that shape the
physics of the ISM of a galaxy, initially studied by McKee & Ostriker (1977)
and followed up by many groups (e.g. Draine, 2011; Kim & Ostriker, 2015; Walch
& Naab, 2015; Gatto et al., 2015; Martizzi et al., 2015; Haid et al., 2016;
Lucas et al., 2020; Bieri et al., 2023; Hirashima et al., 2023). Numerically,
incorporating all the relevant physics in one model for star formation and
feedback has long been a challenge.

The first simulations that include star formation considered a simple over-
density criterion of a single fluid gas, ignoring the multi-phase nature of
the ISM (White & Rees, 1978; Cen & Ostriker, 1992; Navarro & White, 1993; Katz
et al., 1996). Despite the simplicity of these models they have been used to
understand the formation of galaxies in general (e.g., Katz & Gunn, 1991;
Navarro & White, 1994; Steinmetz & Muller, 1995; Mihos & Hernquist, 1996).
However, with the single fluid approach, these models resulted in the well-
known over-cooling problem due to the lack of stellar feedback that resulted
in an overproduction of stars (see, e.g., Somerville & Davé, 2015; Naab &
Ostriker, 2017, for a detailed discussion).

The feedback of massive stars (i.e. stellar explosions, stellar winds, etc.)
is a crucial component of the ISM due to the metal enrichment, builds up
turbulence in the medium and disrupts cold clouds in the neighbour of a star
(see, e.g., Elmegreen & Scalo, 2004; Scalo & Elmegreen, 2004; Naab & Ostriker,
2017, for reviews). Stellar feedback and radiative cooling processes in the
ISM give rise to a multi-phase medium (McKee & Ostriker, 1977) that consists
of hot gas, generated by stellar feedback, warm gas at the cooling/heating
equilibrium and cold gas in the cooling dominated regime that can be split in
neutral and molecular gas, from which the latter is responsible for the
formation of stars. This multi-phase nature of the ISM makes it necessary to
modify the star formation recipes and incorporate the multi-phase structure of
the ISM in the numerical modelling (e.g., Springel & Hernquist, 2003).

However, in addition to the thermal component of the ISM, other factors
influence its structure, such as turbulence, cosmic rays, and magnetic fields.
These ingredients are fundamental components of the energy budget of the ISM
but are very often ignored, especially when it comes to large-scale structure
formation simulations or are not coupled to the sub-grid models for star
formation or ISM modelling.

Recently, a number of different star formation recipes have been developed to
capture processes that are relevant to the multi-phase structure of the ISM.
These models can efficiently be coupled to the heating from supernovae (SNe)
and account for a treatment of non-equilibrium cooling (e.g., Krumholz &
McKee, 2005; Hopkins et al., 2012; Hu et al., 2017; Ebagezio et al., 2023).
Furthermore, there is a class of models that incorporate turbulence (e.g.,
Federrath & Klessen, 2012; Semenov et al., 2018; Kretschmer & Teyssier, 2020)
and models that shift the modelling from a density threshold that initiates
star formation to a pressure threshold (e.g., Murante et al., 2010). Although
the models mentioned above present a significant advancement in the treatment
of the ISM, they still miss more complicated ingredients, such as the impact
of magnetic fields and cosmic rays, which can be easily accounted for in
pressure-based models, (but see, e.g., Girma & Teyssier, 2024, for an
extension of the turbulence models to account for magnetic fields).

Magnetic fields are present on every scale of the Universe and are significant
when ionised matter is involved. On sub-galactic scales, magnetic fields are
the main reason for cosmic ray acceleration and propagation via magnetised
shocks and are relevant for star formation (e.g., Beck, 2015, for a review).
Recent studies (e.g., Pakmor & Springel, 2013; Steinwandel et al., 2020) show
that they could potentially contribute to the formation of galactic outflows.
Observations of galaxies (e.g., Beck & Wielebinski, 2013; Heesen et al., 2023)
indicate the existence of an ordered magnetic field structure with a magnitude
of a few $\mathrm{\mu}$Gauss and a turbulent field of similar strength. It has
also been confirmed that the magnetic energy density of galaxies is of the
same order as the thermal component of the ISM (e.g., Tabatabaei et al., 2008;
Basu & Roy, 2013; Manna & Roy, 2023a). The ISM is a low-$\beta$ plasma (where
$\beta$ is the ratio of the thermal pressure over the magnetic pressure), and
therefore, the magnetic fields can be dynamically important for the global
evolution of galaxies, and it is fundamental that they are accurately modelled
in numerical simulations.

Here, we present an updated version of the widely used ISM model for
cosmological simulations by Springel & Hernquist (2003), taking into account
the effect of the magnetic fields and allowing for the effect of other non-
thermal components of the ISM, such as cosmic rays. This model is particularly
designed for the next generation of cosmological simulations and systems that
usually do not resolve the details of the ISM of galaxies.

This work aims to include the non-thermal effects of the ISM into the model
for a self-consistent treatment of star formation in the field of galaxy
formation and evolution. In particular, our updated model includes the
molecular fraction of the multi-phase ISM, which is responsible for the star-
forming gas, in addition to the cold and hot gas that was included in Springel
& Hernquist (2003). We use observational-driven relations between the
molecular fraction of the gas and the total pressure of the galaxy (Blitz &
Rosolowsky, 2006), which allows us to include the effect of the magnetic
fields via the magnetic pressure as an additional regulation for star
formation. In future studies, other non-thermal pressure components can be
included with our approach. Apart from including magnetic field effects in the
star formation process by taking into account the magnetic pressure in the
formation of cold clouds (i.e., molecular fraction of the ISM), we incorporate
the effect of the stellar feedback on the magnetic field seeding of the
galaxy. This model is proposed in Beck et al. (2013) and assumes that SNe
explosions seed the magnetic field of a galaxy with no need for a pre-existing
primordial field. We couple our novel star-formation model with the one
outlined in Beck et al. (2013) to accurately capture the mutual influence
between star-formation processes and magnetic fields in a self-consistent
manner.

The paper is structured as follows. Sec. 2 describes the equations of the
updated star formation model and the motivation for the selection of free
parameters. Sec. 3 briefly describes the numerical methods and the simulation
setup for an isolated Milky Way-type galaxy simulation embedded in a realistic
environment to test our updated star formation and feedback implementation. In
Sec. 4, we show the results of the Milky Way-type galaxy simulation and
compare the properties of the magnetic field-star formation rate correlation
to observations. We summarize and give conclusions in Sec. 5.

##  2 Multiphase model for star formation and feedback

We present an updated model for star formation and feedback that addresses the
multi-phase structure of the ISM and incorporates magnetic fields to regulate
the star formation process. Our model is an extension of the widely-used star
formation model by Springel & Hernquist (2003) (SH03 hereafter). In addition,
we couple our model to the magnetic seeding model of Beck et al. (2013), which
describes the magnetic field seeding through SNe explosions, which rate is
derived from our updated star formation recipe. Therefore, we self-
consistently model the effect of the magnetic fields in the star formation
process of the galaxy that shapes the ISM and also the effect of the star
formation on the magnetic field evolution in galaxies. In section 2.1, we
present the equation of our model; in section 2.2, we specify and motivate the
selection of the free parameters in our model, and in section 2.3, we
summarize the work of Beck et al. (2013) and how it is altered based on our
model details.

###  2.1 Equations of the model

The new star formation model presented here is a modification of the widely
used SH03 model for star formation and feedback and is implemented in the
developer’s version of the smoothed particle magneto-hydrodynamics (SPMHD)
code gadget (Springel et al., 2001; Springel, 2005; Dolag & Stasyszyn, 2009).
In our model, a resolution element representing a segment of the ISM
encompasses both hot and cold gas. A proportion of the cold gas undergoes a
transition to molecular form, which subsequently gives rise to star formation.
A unique aspect of our model is incorporating the molecular fraction, which is
regulated by pressure, including thermal and non-thermal components of the
ISM. The pressure-regulated molecular fraction of the ISM is absent in SH03.
Therefore, we modify their model accordingly and derive the equations that
describe the sub-grid model for star formation and feedback in this section.

Our model includes thermal heating of SNe that evaporates the surrounding cold
clouds, forms the hot volume filling phase of the ISM and enriches the ambient
ISM with metals. In addition, the hot phase is levelled down due to radiative
cooling, quickly restoring the cold phase that eventually transfers into the
star-forming molecular phase, often neglected in sub-grid modelling for star
formation, which makes the modelling more complex compared to previous
approaches. Once the gas reaches a critical density threshold, it will enter
the star-forming regime. Even though our model still relies on a density
threshold for the gas to undergo star formation, this density threshold is
regulated by the total pressure of the galaxy that includes thermal and non-
thermal components of the ISM, as it will be shown later in this section.

A schematic description of the new model can be found in figure 1. The main
difference compared to the model of SH03 is altering the multi-phase model in
a way that includes an additional component of the ISM, that one of the
molecular form. The molecular form of the ISM gas is directly correlated to
the total pressure, which involves thermal and non-thermal components of the
ISM pressure. We use an observationally driven relation between the pressure
and the molecular fraction of a galaxy by Blitz & Rosolowsky (2006), which
will be further motivated later in this section.

Figure 1: Schematic representation of the mass flow (of a resolution element)
between different phases of the ISM and the interaction between them. Hot gas
is cooling into the cold phase, from where it transits into the molecular
regime before it is turned into stars once it reaches a density threshold. The
stars explode as SNe, which restores the hot ISM phase. Furthermore, the hot
phase interacts with the cold phase as cold clouds evaporate. The novelty of
this model is to include the molecular phase of the ISM, regulated by thermal
and non-thermal pressure components of the ISM, which gives rise to star
formation.

We follow the modelling of SH03 and present here the updated version of the
model equations in our three-phase ISM model. In the updated model, the stars
are formed on a characteristic timescale that we choose to be the free fall
timescale of the cold gas with a star formation efficiency of $f_{*}$.
Therefore, the star formation rate density is described by the following
equation

$\dfrac{d\rho_{*}}{dt}=(1-\beta)f_{*}f_{\mathrm{mol}}\dfrac{\rho_{c}}{t_{\mathrm{dyn}}},$
(1)

with $\rho_{*}$ being the density of the star-forming gas, $f_{\mathrm{mol}}$
is the fraction of the gas in molecular form that is star-forming, $\rho_{c}$
111In all equations of the model, a subscript c refers to the quantities
describing the cold gas of the medium and h refers to the hot gas. is the
density of the cold gas, $t_{\mathrm{dyn}}$ is the dynamical timescale, and
$\beta$ is the fraction of the massive, short-lived stars, explode as SN type
II in timescales comparable to the dynamical time of the system and, therefore
cannot be accounted for the star formation. The parameter $\beta$ depends on
the initial mass function (IMF), and here we choose $\beta=0.1$ for a common
Salpeter IMF (Salpeter, 1955). The star formation rate relates to the gas
reservoir of the galaxy through the observational Schmidt-Kennicutt relation
(Schmidt, 1959; Kennicutt, 1998):

$\Sigma_{\mathrm{SFR}}=(2.5\pm 0.7)\cdot
10^{-4}\left(\dfrac{\Sigma_{\mathrm{gas}}}{\mathrm{M_{\odot}pc^{-2}}}\right)^{1.4\pm
0.15}\mathrm{\dfrac{M_{\odot}}{yr\cdot kpc^{2}}}.$ (2)

We assume that star formation is associated with the dynamical time of the
cold gas, which is the free fall time of the cold clouds given by

$t_{\mathrm{dyn}}=\sqrt{\dfrac{3\pi}{32G\rho_{c}}}.$ (3)

Hereby, $G$ denotes the gravitational constant, and $\rho_{\mathrm{c}}$ is the
(local) density of the cold clouds.

The star formation rate is proportional to the fraction of the molecular
clouds $f_{\mathrm{mol}}$. The molecular fraction of the gas is calculated
given the observational relation found in the study of Blitz & Rosolowsky
(2006). The molecular fraction is proportional to the hydrostatic pressure of
the galaxy and is given by

$\displaystyle f_{\mathrm{mol}}=\dfrac{1}{1+P_{0}/P_{\mathrm{ext}}},$ (4)

with $P_{0}/k_{B}=35000\,\mathrm{Kcm^{-3}}$ and expresses the external
pressure of the ISM when half of the gas is in molecular form and where
$k_{B}$ is the Boltzmann constant. The external pressure $P_{\mathrm{ext}}$ is
calculated from the mid-plane pressure of an infinite disk that consists of
gas and stars, assuming that the gas scale height is much lower than the
stellar scale height. The pressure is derived from the following expression as
shown in Blitz & Rosolowsky (2006)

$P_{\mathrm{ext}}=(2G)^{0.5}\Sigma_{g}u_{g}\left[\rho_{\mathrm{stars}}^{0.5}+\left(\dfrac{\pi}{4}\rho_{g}\right)\right].$
(5)

In the last equation, $\Sigma_{g}$ is the surface density of the gas, $u_{g}$
is the vertical velocity dispersion of the gas, $\rho_{\mathrm{stars}}$ is the
mid-plane density of the stellar population and $\rho_{g}$ the mid-plane
density of the gas. This pressure corresponds to the (total) hydrostatic
pressure of the galaxy, which, apart from the thermal pressure, also includes
non-thermal components, such as the magnetic field pressure. For the
simulations in this work, we approximate this pressure with the SPMHD pressure
of the particle plus the magnetic pressure given by the SPMHD code gadget.

From the stars that are formed, many of them explode as SNe, which heats the
surrounding medium in the form of thermal feedback. Thus, the hot gas gains
energy at the rate

$\left.\dfrac{d}{dt}(\rho_{h}u_{h})\right|_{\mathrm{SN}}=\epsilon_{\mathrm{SN}}f_{*}f_{\mathrm{mol}}\dfrac{d\rho_{*}}{dt}=\beta
u_{\mathrm{SN}}f_{*}f_{\mathrm{mol}}\dfrac{d\rho_{c}}{t_{\mathrm{dyn}}},$ (6)

with $\epsilon_{\mathrm{SN}}$ the energy per unit mass that corresponds to one
SN explosion $\epsilon_{\mathrm{SN}}=4\cdot 10^{48}\,\mathrm{erg\cdot
M_{\odot}^{-1}}$ for the IMF adopted in this description (Salpeter, 1955). In
the last equation, $\rho_{h}$ and $u_{h}$ refer to the density and energy per
unit mass of the hot phase, respectively.

Figure 2: Fraction of the cold phase of the ISM ($x$ in the model equations)
as a function of density. The orange points show the result of the new
pressure-based star formation model, and the green points represent the
outcome of the classic Springel & Hernquist (2003) model for star formation.

As a result of the heating, the SN blast wave destroys the star-forming clouds
around the region of the SN explosion. Thus, the cold phase evaporated, and
the cold mass fraction of the ISM loses mass at the rate

$\left.\dfrac{d\rho_{c}}{dt}\right|_{\mathrm{EV}}=A\beta
f_{*}f_{\mathrm{mol}}\dfrac{\rho_{c}}{t_{\mathrm{dyn}}},$ (7)

with $A$ being the evaporation parameter that scales with density as
$A\sim\rho^{-4/5}$ (McKee & Ostriker (1977), see also section 2.2). The growth
of the cold clouds comes from the radiative cooling of the hot gas. We follow
the formulation of SH03 that assumes the thermal instability (TI) operating,
giving rise to cold clouds due to the cooling of hot gas. The mass flow from
one phase to the other due to radiative cooling is described as

$\left.\dfrac{d\rho_{c}}{dt}\right|_{\mathrm{TI}}=-\left.\dfrac{d\rho_{h}}{dt}\right|_{\mathrm{TI}}=\dfrac{1}{u_{h}-u_{c}}\Lambda_{\mathrm{net}}(\rho_{h},u_{h}),$
(8)

with $u_{c}$ being the energy per unit mass of the cold phase, $u_{h}$ the
energy per unit mass of the hot phase and
$\Lambda_{\mathrm{net}}(\rho_{h},u_{h})$ the cooling function which is
computed from radiative processes taking place in a primordial plasma of
Hydrogen and Helium, as presented by Katz et al. (1996). In gadget, we include
collisional excitation of H and He+, ionization of H, He, He+, recombination
of H+, He+, He++ and free-free (thermal Bremsstrahlung) emission. The gas
cannot cool below $\simeq 10^{4}\,\mathrm{K}$ because lower temperatures would
require proper treatment of molecular cooling, which can be included in our
model but is currently neglected. Thus, we set the temperature of cold clouds
to remain constant at $T_{c}=10^{3}\,\mathrm{K}$.

Taking into consideration all the processes being described, the density of
the hot and cold phases changes as

$\displaystyle\dfrac{d\rho_{c}}{dt}$
$\displaystyle=-f_{*}f_{\mathrm{mol}}\dfrac{\rho_{c}}{t_{\mathrm{dyn}}}-A\beta
f_{*}f_{\mathrm{mol}}\dfrac{\rho_{c}}{t_{\mathrm{dyn}}}+\dfrac{1-f}{u_{h}-u_{c}}\Lambda_{\mathrm{net}}(\rho_{h},u_{h}),$
(9) $\displaystyle\dfrac{d\rho_{h}}{dt}$ $\displaystyle=\beta
f_{*}f_{\mathrm{mol}}\dfrac{\rho_{c}}{t_{\mathrm{dyn}}}+A\beta
f_{*}f_{\mathrm{mol}}\dfrac{\rho_{c}}{t_{\mathrm{dyn}}}-\dfrac{1-f}{u_{h}-u_{c}}\Lambda_{\mathrm{net}}(\rho_{h},u_{h}).$
(10)

In both differential equations, the first term on the right-hand side
represents the gain/loss from the star formation, the second term accounts for
the cold cloud evaporation, and the last one is the effect of the thermal
instability. The parameter $f$, which can have a value of 1 or 0, represents
the onset of the thermal instability. For $f=0$ the thermal instability is
operating, and stars are forming. In the opposite case, when $f=1$, ordinary
cooling takes place. Following SH03, to physically differentiate the two
cases, we use a density threshold, $\rho_{\mathrm{thr}}$, which will be
calculated in section 2.2. Thus, star formation takes place for regions where
the density of the gas exceeds the density threshold $\rho_{\mathrm{thr}}$,
i.e. when $\rho>\rho_{\mathrm{thr}}$.

The energy of the total gas changes according to all the processes mentioned
above and is described in the equation

$\begin{split}\dfrac{d(\rho_{h}u_{h}+\rho_{c}u_{c})}{dt}=-\Lambda_{\mathrm{net}}(\rho_{h},u_{h})+&\beta
f_{*}f_{\mathrm{mol}}\dfrac{\rho_{c}}{t_{\mathrm{dyn}}}u_{\mathrm{SN}}-\\\
&(1-\beta)f_{*}f_{\mathrm{mol}}\dfrac{\rho_{c}}{t_{\mathrm{dyn}}}u_{c}.\end{split}$
(11)

As mentioned before, the temperature of the cold clouds remains constant at
$T_{c}=10^{3}\,\mathrm{K}$. Consequently, we assume that $u_{c}$ is constant
in the above equations, and therefore, we can follow the evolution of the
specific energy of the hot phase.

The hot phase of the gas will evolve according to the differential equation

$\begin{split}\dfrac{du_{h}}{dt}&=-\left(\beta\dfrac{f_{*}f_{\mathrm{mol}}}{t_{\mathrm{dyn}}}\dfrac{\rho_{c}}{\rho_{h}}+\dfrac{A\beta}{t_{\mathrm{dyn}}}\dfrac{\rho_{c}}{\rho_{h}}\right)u_{h}+\\\
&f_{*}f_{\mathrm{mol}}\dfrac{\beta}{t_{\mathrm{dyn}}}\dfrac{\rho_{c}}{\rho_{h}}(u_{\mathrm{SN}}+u_{c})+f_{*}f_{\mathrm{mol}}\dfrac{A\beta}{t_{\mathrm{dyn}}}\dfrac{\rho_{c}}{\rho_{h}}u_{c}.\end{split}$
(12)

The last equation results in an equilibrium solution for timescales larger
than the dynamical time. Thus, the temperature of the hot phase will evolve
towards an equilibrium state, as described by

$u_{h}=\dfrac{u_{\mathrm{SN}}}{A+1}+u_{c}.$ (13)

Deviations from this equilibrium solution are decaying on a timescale

$\tau_{h}=\dfrac{t_{\mathrm{dyn}}\rho_{h}}{\beta(A+1)f_{*}f_{\mathrm{mol}}\rho_{c}}.$
(14)

It is interesting to notice that this timescale depends not only on the cold
fraction of the gas but also on the molecular fraction of the cold gas.
Therefore, when the magnetic fields are amplified, the molecular fraction
rises (as the total pressure rises) and adjusts the equilibrium state compared
to the one that does not include magnetic fields. The growth of cold clouds is
balanced by the star formation and SN feedback towards an equilibrium which
self-regulates these processes during galaxy evolution. The gas then behaves
as an effective medium with pressure

$P_{\mathrm{eff}}=(\gamma-1)(\rho_{h}u_{h}+\rho_{c}u_{c}).$ (15)

The effective pressure remains constant due to equilibrium. Taking the above
into account equation 11 yields

$\dfrac{\rho_{c}}{t_{\mathrm{dyn}}}=\dfrac{\Lambda_{\mathrm{net}}(\rho_{h},u_{h})}{f_{*}f_{\mathrm{mol}}(\beta
u_{\mathrm{SN}}-(1-\beta)u_{c})}.$ (16)

We can derive the cold fraction of the gas as a function of the gas density.
The cold fraction is defined as $x=\dfrac{\rho_{c}}{\rho}$. The cooling
function of the hot phase is
$\Lambda_{\mathrm{net}}(\rho_{h},u_{h})=(\rho_{h}/\rho)^{2}\Lambda_{\mathrm{net}}(\rho,u_{h})$.
Using the definition

$y(x)=\dfrac{t_{\mathrm{dyn}}(x)\Lambda_{\mathrm{net}}(\rho_{h},u_{h})}{f_{*}f_{\mathrm{mol}}(\beta
u_{\mathrm{SN}}-(1-\beta)u_{c})},$ (17)

the ratio of the densities of the hot and cold phases is written as

$\dfrac{\rho_{c}}{\rho_{h}}=\dfrac{\rho_{h}}{\rho}y(x).$ (18)

As a result, the cold fraction can be written as a function of the parameter
$y$ as

$x=1+\dfrac{1}{2y(x)}-\sqrt{\dfrac{1}{y(x)}+\dfrac{1}{4y(x)^{2}}}.$ (19)

In our model, the cold fraction does not depend only on $y$, but on the cold
fraction itself since the parameter $y$ is a non-linear function of $x$. The
non-linear expression for the cold fraction comes from choosing the dynamical
time of the cold gas as a characteristic timescale of the star formation. The
dynamical time (see equation 3) depends on the cold fraction $x$ as

$t_{\mathrm{dyn}}=x^{-1/2}\sqrt{\dfrac{3\pi}{32G\rho}}.$ (20)

Given the non-linearity of the dependencies, we can not solve for $x$ but
rather find its value numerically. The cold fraction is the root of the
function

$f(x)=1+\dfrac{\sqrt{x}}{2\delta}-\sqrt{\dfrac{x}{\delta}+\dfrac{x}{2\delta^{2}}}-x,$
(21)

with $\delta$

$\delta=\sqrt{\dfrac{3\pi}{32G\rho}}\cdot\dfrac{\Lambda_{\mathrm{net}}(\rho_{h},u_{h})}{f_{*}f_{\mathrm{mol}}\rho(\beta
u_{\mathrm{SN}}-(1-\beta)u_{c})}.$ (22)

After testing that the equation 21 has a solution, we numerically solve it
with a simple bisection method. The dependence of the cold fraction $x$ on the
gas density is shown in figure 2.

For comparison with the star formation model of SH03, we directly compare our
updated version to the original model in the same figure. The obvious
difference is that, in our case, the cold gas fraction reaches lower values
and also does not converge to 1 for high densities. This comes from the
dependence of the dynamical time on the cold fraction that inserts a
complicated correlation between the cold gas fraction $x$ and the total gas
density $\rho$. Thus, for high densities, we don’t allow all the gas to be in
the cold phase and form more stars, but we have included other regulation
factors such as the external pressure of the ISM through the molecular
fraction $f_{\mathrm{mol}}$.

In order to distinguish between the different states of the model, i.e. star-
forming or non-star-forming, the gas parcel enters the multi-phase model (i.e.
star-forming) once its density is higher than a density threshold,
$\rho_{\mathrm{thr}}$ (see discussion on equation 9 and equation 10). Below
this density, we set the cold fraction to $x=0$, which means that the gas is
not star-forming. The density threshold is calculated from the model equations
and will be discussed in the next section.

###  2.2 Selection of Parameters

Figure 3: Star formation rate surface density as a function of the surface
density of the total gas in a galaxy for different star formation efficiency
parameters. The black solid line is the observational Schmidt-Kennicutt
relation, and the grey lines represent its errors. All quantities are
azimuthically averaged.

The evaporation factor $A$ (introduced in equation 7) scales with the local
density of the ISM, which is theoretically motivated by McKee & Ostriker
(1977). We scale it with the density threshold as

$A(\rho)=A_{0}\left(\dfrac{\rho}{\rho_{\mathrm{thr}}}\right)^{-4/5}.$ (23)

In order to constrain $A_{0}$, we follow SH03 and consider the onset of the
thermal instability. From the equilibrium solution in equation 13, the
temperature of the hot phase is approximately $10^{5}$ K, which is the peak of
the cooling function; therefore, the cooling is very effective, and the
thermal instability is triggered. Thus $T_{\mathrm{SN}}/A_{0}=10^{5}$ K, which
constrains the parameter $A_{0}$ to be $A_{0}=1000$, for a typical SN
temperature of $T_{\mathrm{SN}}=2\mu
u_{\mathrm{SN}}m_{\mathrm{H}}/3k_{\mathrm{B}}=10^{8}$ K.

In the following, we derive the density threshold, which differentiates
between star-forming and non-star-forming gas using equation 17. The cooling
function can be written as
$\Lambda_{\mathrm{net}}(\rho,u_{h})/\rho=u_{h}/t_{\mathrm{cool}}$. In order to
calculate the density threshold, we calculate the value of the cooling
function at the baryon over-density calculated as
$\rho_{\mathrm{ov}}=10^{5}\dfrac{3H^{2}}{8\pi G}$, following SH03. Therefore,
the density threshold is written as

$\rho_{\mathrm{thr}}=\dfrac{x_{\mathrm{thr}}}{(1-x_{\mathrm{thr}})^{2}}\dfrac{f_{\mathrm{mol}}(\beta
u_{\mathrm{SN}}-(1-\beta)u_{c})}{f_{*}t_{\mathrm{dyn}}(\rho_{\mathrm{ov}})\Lambda(\rho_{\mathrm{ov}},u_{\mathrm{SN}}/A_{0})}.$
(24)

The cooling function is calculated at the threshold where
$u_{h}=u_{\mathrm{SN}}/A_{0}+u_{c}\simeq u_{\mathrm{SN}}/A_{0}$. The cold
fraction at the threshold $x=x_{\mathrm{thr}}$ is given by
$x_{\mathrm{thr}}=1+(A_{0}+1)(u_{c}-u_{4})/u_{\mathrm{SN}}\simeq
1-A_{0}u_{4}/u_{\mathrm{SN}}$. The latter quantity is calculated by setting
the condition $u_{\mathrm{eff}}(\rho_{\mathrm{thr}})=u_{4}$, where $u_{4}$ is
the specific energy that corresponds to temperature of
$T_{4}=10^{4}\,\mathrm{K}$. This implies that the pressure is a continuous
function of density and should stay constant before and after the onset of the
star formation, which occurs at $\rho=\rho_{\mathrm{thr}}$. A further step
would require considering the magnetic pressure of other non-thermal pressure
components at the threshold. To be more specific, instead of taking into
account only the effective specific energy of the gas, we could include the
magnetic energy or alternatively consider the sum of the magnetic pressure
plus the gas pressure being a continuous function of density. The magnetic
pressure and energy scales with $B^{2}$, which scales with the density as
$\rho^{4/3}$ in the flux freezing limit. This would include more complications
in calculating the cold fraction at the threshold and is not crucial for the
simple approach we are adopting here.

In equation 24, we notice that the density threshold depends on the molecular
fraction of the gas. Since we have already set the condition of continuous
pressure on the onset of star formation to be $P_{4}$, i.e., the pressure that
corresponds to a temperature of $10^{4}\,\mathrm{K}$, it is then a natural
consequence to use this pressure to calculate the molecular fraction at the
threshold, which will be

$f_{\mathrm{mol}}(P_{4})=\dfrac{1}{1+P_{0}/P_{4}}.$ (25)

The density threshold is now calculated as

$\rho_{\mathrm{thr}}=\dfrac{x_{\mathrm{thr}}}{(1-x_{\mathrm{thr}})^{2}}\dfrac{f_{\mathrm{mol}}(P_{4})(\beta
u_{\mathrm{SN}}-(1-\beta)u_{c})}{f_{*}t_{\mathrm{dyn}}(\rho_{\mathrm{ov}})\Lambda(\rho_{\mathrm{ov}},u_{\mathrm{SN}}/A_{0})}.$
(26)

Another parameter that we should specify, which is a free parameter of the
model, is the efficiency of the star formation process $f_{*}$. The efficiency
of star formation is related to the gas component of the galaxy, especially
the cold gas. These two quantities are correlated, as observations show,
according to the Schmidt-Kennicutt relation (Schmidt, 1959; Kennicutt, 1998,
see equation 2), which is a tight relation between the star formation of the
galaxy and the surface density of the gas in the galaxy. In order to constrain
the efficiency parameter $f_{*}$ we run test simulations of an isolated Milky
Way type galaxy, with different $f_{*}$ as shown in figure 3. In order to
reproduce the observed Schmidt-Kennicutt relation, the value of the efficiency
parameter should be $f_{*}=0.1$, i.e. star formation efficiency of $10\%$.

Moreover, as far as the molecular fraction is concerned, we set the parameter
$P_{0}$ to be $P_{0}/k_{\mathrm{B}}=35000\,\mathrm{Kcm^{-3}}$ according to the
observations presented in Blitz & Rosolowsky (2006). The authors distinguish
two different groups of galaxies, and the parameter $P_{0}$ changes
significantly for each of them. There are three galaxies that on average they
find $P_{0}/k_{\mathrm{B}}=7700\,\mathrm{Kcm^{-3}}$ and for the rest of the
galaxies they find $P_{0}/k_{\mathrm{B}}=43000\,\mathrm{Kcm^{-3}}$. They
attribute this difference to the content of neutral hydrogen that each galaxy
contains. As calculated from observations, the mean value of the projected
surface density of the neutral hydrogen may be low due to intense tidal or ram
pressure stripping. However, when testing our model with the different values
of this parameter the one that follows the Schmidt-Kennicutt relation remains
the $P_{0}$/$k_{\mathrm{B}}=35000\,\mathrm{Kcm^{-3}}$.

The last parameter that we examine is the temperature of cold clouds. This was
assumed to be constant at $T_{c}=1000\,\mathrm{K}$, but this is not close to
the temperatures at which stars are formed in nature. Therefore, we examine
the case of lowering the cold gas temperature and run a test simulation for
$T_{c}=300\,\mathrm{K}$, which is motivated by Murante et al. (2015). The
lower temperature barely affects the cold fraction and the overall star
formation behaviour. The density threshold changes slightly because of its
dependence on $u_{c}$. However, we should bear in mind that the gas at
temperatures lower than $10^{3}-10^{4}\,\mathrm{K}$ needs molecular cooling
for proper treatment. Since this is not included in our simulations, we will
keep the temperature of the cold clouds constant to $10^{3}\,\mathrm{K}$. A
summary of the model parameters is shown in table 1.

Table 1: Summary of the parameters used in the star formation model, where $f_{*}$ is the star formation efficiency, $P_{0}/k$ expresses the external pres- sure of the ISM, $\beta$ is the fraction of the massive, short-lived stars that explode as core-collapse SN, and $A_{0}$ is the scaling factor for the evaporation factor $A$.  parameter | $f_{*}$ |  $P_{0}/k$ $[\mathrm{Kcm^{-3}}]$ | $\beta$ | $A_{0}$ |  $T_{c}$ $[\mathrm{K}]$
---|---|---|---|---|---
value | 0.1 | 35000 | 0.1 | 1000 | 1000

###  2.3 Supernova Magnetic field seeding

For a self-consistent model, we also aim to include the effects of star
formation via the stellar feedback on the magnetic field evolution of
galaxies. To do so, we couple the star formation model with the SN magnetic
seeding model of Beck et al. (2013). This model assumes that the magnetic
fields do not have a primordial origin but are deposited in the ISM from SN
explosions. Here, we give a summary of the model by Beck et al. (2013) and how
this is altered because we couple it to the updated star formation model
discussed in the previous sections.

According to Beck et al. (2013), an extra term is included in the induction
equation to include the effect of the magnetic field seeding. Thus, the
induction equation is written as

$\dfrac{\partial\mathbf{B}}{\partial
t}=\nabla\times(\mathbf{v}\times\mathbf{B})+\eta\nabla^{2}\mathrm{B}+\left.\dfrac{\partial\mathbf{B}}{\partial
t}\right|_{\mathrm{seed}},$ (27)

with the seeding term

$\left.\dfrac{\partial\mathbf{B}}{\partial
t}\right|_{\mathrm{seed}}=\sqrt{N_{\mathrm{SN}}^{\mathrm{eff}}}\dfrac{B_{\mathrm{inj}}}{\Delta
t}\mathbf{e}_{B},$ (28)

where $\mathbf{e}_{B}$ is a unit vector in the direction of the seeding,
$B_{\mathrm{inj}}$ is the amplitude of the seeding field and
$N_{\mathrm{SN}}^{\mathrm{eff}}$ is a normalization constant that is connected
with the number of SN explosions. The last parameter is not a free parameter
of the model but is calculated directly from the star formation recipe in
gadget. For the updated star formation model, the mass of stars, $m_{*}$, that
are formed in each timestep $\Delta$t is calculated as

$m_{*}=f_{*}f_{\mathrm{mol}}m_{c}\dfrac{\Delta t}{t_{\mathrm{dyn}}}$ (29)

where $m_{c}$ is the mass of the cold clouds. The effective number of SNe
explosions is given by

$N_{\mathrm{SN}}^{\mathrm{eff}}=\alpha m_{*}$ (30)

with $\alpha$ being a parameter that specifies the number of SNe explosions
per solar mass. For our case the parameter $\alpha$ is
$\alpha=0.008\,\mathrm{M_{\odot}^{-1}}$ (for a Salpeter IMF).

The total injected magnetic field for all SN explosions is given by

$B_{\mathrm{inj}}^{\mathrm{all}}=\sqrt{N_{\mathrm{SN}}^{\mathrm{eff}}}B_{\mathrm{SN}}\left(\dfrac{r_{\mathrm{SN}}}{r_{\mathrm{SB}}}\right)^{2}\left(\dfrac{r_{\mathrm{SB}}}{r_{\mathrm{inj}}}\right)^{3}$
(31)

where $B_{\mathrm{SN}}$ is the mean strength of magnetic fields in SN
explosions, $r_{\mathrm{SB}}$ is the radius of the super-bubble of the
explosion assuming spherical geometry for the remnant, $r_{\mathrm{SN}}$ is a
typical radius of the SN remnant and $r_{\mathrm{inj}}$ is determined by the
smoothing length of the simulation and shows the region where the super-
bubbles are placed. The magnetic field seeding rate for this model is
calculated as

$\dot{B}_{\mathrm{seed}}\simeq
B_{\mathrm{SN}}\left(\dfrac{r_{\mathrm{SN}}}{r_{\mathrm{SB}}}\right)^{2}\left(\dfrac{r_{\mathrm{SB}}}{r_{\mathrm{inj}}}\right)^{3}\dfrac{\sqrt{\dot{N_{\mathrm{SN}}}\Delta
t}}{\Delta t}.$ (32)

The SN radius is assumed to be $r_{\mathrm{SN}}=5\,\mathrm{pc}$, a typical
value for the magnetic field strength is $B_{\mathrm{SN}}=10^{-4}\,\mathrm{G}$
which is then distributed in a bubble of radius
$r_{\mathrm{SB}}=25\,\mathrm{pc}$.

The magnetic field configuration added to the induction equation should be
divergence-free. Therefore, a straightforward way is to assume a dipole
structure of the magnetic field seed, so the seeding field changes at a rate

$\left.\dfrac{\partial\mathbf{B}}{\partial
t}\right|_{\mathrm{seed}}=\dfrac{1}{{\mathbf{r}}^{3}}\left[3\left(\dfrac{\partial\mathbf{m}}{\partial
t}\cdot\mathbf{e}_{r}\right)\mathbf{e}_{r}-\dfrac{\partial\mathbf{m}}{\partial
t}\right]$ (33)

with $\mathbf{m}$ the dipole magnetic moment, $\mathbf{e}_{r}$ is a unit
vector in the $\mathbf{r}$ direction. The time derivative of each dipole
moment is written as

$\dfrac{\partial\mathbf{m}}{\partial
t}=\sigma\dfrac{B_{\mathrm{inj}}^{\mathrm{all}}}{\Delta t}\mathbf{e}_{B}.$
(34)

Here $\mathbf{e}_{B}$ is the direction of the seed magnetic field, which is
chosen to be in the direction of the acceleration of the particle, hence
$\mathbf{e}_{B}=a/\left\lvert a\right\rvert$. The normalization constant
$\sigma$ is given by

$\sigma=r_{\mathrm{inj}}^{3}\sqrt{\dfrac{1}{2}\bar{f}^{3}(1+\bar{f}^{3})}$
(35)

with $\bar{f}=r_{\mathrm{soft}}/r_{\mathrm{inj}}$ is the ratio between the
softening and injection length. The parameter $\sigma$ is used to normalize
the energy injected, soften the magnetic dipoles in the centre, and truncate
it in the scale of the injection length.

##  3 Numerical Method

We employ the developer’s version of the Tree Smooth Particle
Magnetohydrodynamics (SPMHD) code gadget, which is based on the well-tested
code gadget2 (Springel, 2005) and incorporates an updated SPH implementation
introduced in Beck et al. (2016). The code follows a Lagrangian treatment of
magnetohydrodynamics and a tree algorithm for Newtonian gravity. The magnetic
field implementation is described in detail in Dolag & Stasyszyn (2009). The
aforementioned model (section 2) for star formation and feedback is coupled to
this version of the code and allows a realistic and improved treatment of the
ISM in large-scale simulations. This includes star formation, the SN feedback
from stars, radiative cooling and magnetic fields.

For the initial conditions, we follow the methods for a galactic model first
described in Hernquist (1993) in the numerical implementation that has been
used in Springel (2005). The galactic model consists of a stellar bulge, a
stellar disk, and a gas disk embedded in a dark matter halo. The dark matter
and bulge follow a mass distribution of a Hernquist profile (Hernquist, 1993).
This system is surrounded by a hot circum-galactic medium (CGM) as first
presented in Moster et al. (2010); here, we use the implementation of
Steinwandel et al. (2019) that incorporates some modifications from Donnert
(2014) concerning the sampling mechanism of the CGM’s density distribution.
This follows a $\betaup$-profile, namely

$\rho(r)=\rho_{0}\left(1+\dfrac{r^{2}}{r_{c}^{2}}\right)^{-3\betaup/2}.$ (36)

With $\betaup=2/3$, following Mastropietro & Burkert (2008). In the last
equation, $\rho_{c}=5\cdot 10^{-26}\,\mathrm{g/cm^{3}}$ and $r_{c}=0.33\,$kpc,
as shown in Steinwandel et al. (2019). We chose this setup in order to
accurately model the smooth accretion of gas towards the galactic disk, and as
we use particles to discretize the induction equation, we need a carrier
medium to properly deal with the vacuum boundary conditions of the magnetic
field. This set-up is already tested and described in detail in Steinwandel et
al. (2019). Although this galactic model lacks a cosmological background, it
is a useful layout in order to examine the effects of the updated star
formation routines in a clean environment. In table 2, a summary of the
parameters for the initial conditions used is shown.

The galaxy has mass resolution of $4800\,\mathrm{M_{\odot}}$ in the gas and
stellar component and $96000\,\mathrm{M_{\odot}}$ in dark matter. The spatial
resolution is limited by the gravitational softening, which is calculated from
the recursive formula

$\centering\epsilon=\epsilon_{\mathrm{old}}\left(\dfrac{m}{m_{\mathrm{old}}}\right)^{1/3}.\@add@centering$
(37)

The reference simulation for the softening that we used to obtain the force
softening of our simulations is a set of
magneticum222http://www.magneticum.org/simulations.html simulations
(Hirschmann et al., 2014) as a benchmark, denoted by the subscript “old” in
equation 37. For the gas particles the gravitational softening it is
$\epsilon_{\mathrm{gas}}=10\,\mathrm{pc/h}$, for the stellar particles it is
$\epsilon_{\mathrm{stars}}=20\,\mathrm{pc/h}$ , and for the dark matter
particles it is $\epsilon_{\mathrm{DM}}=83\,\mathrm{pc/h}$.

Table 2: List of parameters for the initial conditions of the galaxy simulation surrounded by a gas halo. Disk Parameters |  |
---|---|---
Total mass [$10^{10}\ \mathrm{M_{\odot}}$] | $M_{200}$ | 100
Virial radius | $r_{200}$ | 145
Halo concentration | $c$ | 12
Spin parameter | $\lambda$ | 0.033
Disk spin fraction | $j_{d}$ | 0.067
Disk mass fraction | $m_{d}$ | 0.067
Bulge mass fraction | $m_{b}$ | 0.034
Disk scale length [kpc] | $l_{d}$ | 2.1
Disk height [ld] | $z_{0}$ | 0.2
Bulge size [ld] | $l_{b}$ | 0.2
Table 3: Total particle number for each particle type of the galaxy, the gas and dark-matter halo system used as initial conditions for our simulations. Particle Numbers [$10^{6}$] |  |
---|---|---
Gas in the disk | $N_{\mathrm{gd}}$ | 1.2
Stellar Disk | $N_{\mathrm{sd}}$ | 4.8
Stellar Bulge | $N_{b}$ | 2.0
Dark matter | $N_{\mathrm{DM}}$ | 6.9

##  4 Results for a Milky Way type galaxy

In this section, we present the results of a Milky Way-type galaxy (High Mass
Galaxy, HMG hereafter), considering two different configurations for the
initial conditions of the magnetic field. In the first model, we use a
primordial magnetic field of strength $10^{-9}\,\mathrm{B}$ in $\hat{x}$
direction of a Cartesian coordinate system (to be called HMG-B), which is
chosen to resemble an already partially amplified magnetic field. The second
model includes the full version of the new star formation model coupled with
the SN seeding prescription (to be called HMG-snB). The magnetic seeding from
SN explosions is a self-consistent way to account for the magnetic fields in
the galaxy and allow the interaction with the star formation and feedback
processes within the ISM. Additionally, we include a run without magnetic
fields, i.e., HMG-noB for comparison. To directly compare with SH03, we run
two different simulations using the SH03 model, i.e., one without magnetic
fields, to be referred to as HMG-noB-SH and one with the full SH03 star
formation and SN feedback model, to be called HMG-snB-SH.

###  4.1 General properties and morphology of the galaxy

Before taking a closer look at the morphology of the galaxy and the structure
of the magnetic field in the galaxy, it is interesting to examine the
temperature-density phase diagram. In Figure figure 4, we show the temperature
of the gas as a function of the density at the start of the simulation for the
model HMG-B for two different points in time, i.e. at the initialisation of
the galaxy simulation (upper panel) and after 1 Gyr of evolution (bottom
panel). The colour bar shows the star formation rate for each computational
particle. The prominent upper-left purple patch corresponds to the hot gas
halo surrounding the galaxy, and it is evident that a part of it cools down as
it accretes onto the galaxy. The bottom-right part of the plot corresponds to
the gas disk of the galaxy. There is a clear separation based on the star
formation rate, which is correlated to the density threshold imposed by our
modified star formation recipe. It is worth mentioning that this density
threshold is correlated with the molecular fraction we obtain directly from
the total hydrostatic pressure in the galaxy’s mid-plane. Moreover, we note
that the temperature is interpreted as an effective temperature of the
particles in the multi-phase model since the cold clouds are kept to a
constant temperature of $T_{c}=10^{3}$ K. This is usually much lower than the
temperature of the hot gas in the ISM, which easily exceeds $10^{5}$ K (McKee
& Ostriker, 1977).

As the galaxy interacts with the CGM, the gas halo cools down, and there is an
inflow of fresh gas onto the galactic disk. A fraction of the gas is not star-
forming with higher temperatures that stream from the cold disc towards the
hot gas halo. This fraction of the gas is pushed out of the disc by the rising
magnetic pressure, which generates an outflow of low, non-star-forming flow of
gas that rises above the disc until it reaches pressure equilibrium with the
surrounding hot gas of the CGM. The gas pushed out remains always below the
density threshold for star formation. We find a very similar structure in the
phase diagrams for the model HMG-snB.

Figure 4: Temperature-density phase diagram, colour-coded by the star
formation rate of each gas particle. The left upper part of the diagram
corresponds to the gas halo particles, while the lower right part corresponds
to the gas content of the galaxy. The top panel is at the beginning of the
simulation, and the bottom is after 1 Gyr of evolution. Figure 5: Projection
on the xy plane of the galaxy for the HMG-B simulation, colour-coded by the
star formation rate. The star formation rate is integrated along the z-axis
(perpendicular to the plane shown). The left panel corresponds to 1 Gyr of
evolution and the right to 2 Gyr. A very similar structure is also found in
the HMG-snB model. Figure 6: Radial profile of the star formation rate of the
galaxy after 1, 2 and 3 Gyr of evolution. The blue lines show the HMG-snB
simulation, and the red lines show the HMG-B model. The prominent peak at 3
Gyr of evolution is due to the outflow created; see text for details.

In figure 5, we show a projection of the star formation rate of the galaxy for
the model HMG-B after 1 Gyr (left panel) and 2 Gyr (right panel) of the
evolution of the system. We find similar morphological structures for the
model HMG-snB. This can be explained by the fact that, at this time, the
magnetic fields are low in both models and, therefore, are not significant
enough to alter the dynamics of the galaxy. The star formation rate follows
the density structure of the galaxy. Thus, our new model predicts higher star
formation rates in the centre and the spiral arms, compared to the inter-arm
regions and the galaxy outskirts, where the star formation rate is generally
lower by a factor of 10. At the beginning of the simulation, the star
formation is prominent in the spiral arms but also in the inter-arm regions.
However, at later times, the star formation concentrates in the spiral arms as
they are highly compressed, and the gas forms stars very efficiently over a
few dynamical times. We also notice that the star formation region shrinks
after 2 Gyr of evolution (see figure 5). This can be interpreted by the fact
that the star formation continuously drops and drops faster than the gas
accretion from the CGM. At later times, star formation is concentrated in the
centre of the galaxy in a bar-like structure fed by the hot gas from the CGM.
As the CGM consists of slowly rotating low angular momentum gas, the hot
material quickly assembles in the centre of the disc, triggering the formation
of a bar that makes the accretion of hot gas to the centre even more
efficient. This process subsequently triggers an outflow that is generated by
the magnetic pressure. We will discuss the process in more detail in section
4.2.

Figure 7: Radial pressure profile for 1,2 and 3 Gyr of evolution. The blue
lines show the gas pressure and the magenta lines show the magnetic field
pressure. For both cases, the solid lines are the pressure components from the
HMG-B simulation, while the dashed lines are from the HMG-snB. Figure 8: Total
star formation rate of the galaxy as a function of time. The blue line shows
the rate for the HMG-snB model and the red line for the HMG-B. The vertical
lines correspond at the start of the outflow for each model following the same
colour code.

###  4.2 Radial evolution

To gauge the consistency of our modified star formation approach, we discuss
radial profiles of the fundamental quantities directly affected by our new
implementation.

In figure 6, the star formation rate over a radius of the galaxy is shown for
the two different models, i.e.HMG-B, HMG-snB, at three different points in
time, i.e., 1 Gyr, 2 Gyr and 3 Gyr. The latter refers to the final time of the
simulation. At early times (t = 1 Gyr), the star formation peaks in the centre
and declines in a smooth power-law as a function of the radius. The
exponential density profile of the galactic disk can easily explain this
behaviour. Both models agree very well and do not lead to substantial
differences in the star formation rate across the disk at early times. The
sharp drop at around 12 kpc occurs due to the boundary between the galactic
disk and the ambient CGM. As low angular momentum gas is accreted towards the
outskirts of the disk, it quickly cools down and forms stars in the outer
parts before the star formation rate drops to zero within the CGM.

At intermediate times (t = 2 Gyr, middle panel), we notice that the star
formation rate is slightly higher in the central part of the galaxy for the
HMG-B model. In that model, the magnetic field pre-exists and thus, at 2 Gyr,
it has already been amplified enough to assist star formation, as suggested by
our model. In contrast, the HMG-snB model starts with no magnetic field, and
therefore, it needs more time to build up a field with sufficient strength.
The same reason applies to the outskirts of the galaxy. In the HMG-snB model,
the star formation rate drops off in a smaller radius compared to the HMG-B,
because there is no magnetic field at initialization, leading to no
contribution of the magnetic field pressure in the molecular fraction of the
gas and therefore in the star formation rate.

The right panel shows the evolution after 3 Gyr. The star formation rate shows
multiple peaks across the galaxy because prominent spiral arms are formed, as
also shown in figure 5. The star formation rate for the HMG-snB model shows
different and more prominent spikes (compared to the HMG-B model) because the
star formation directly follows the magnetic field seeding process in the
spiral arms, which increases the molecular fraction in those regions and
subsequently that increases the star formation rate.

The radial pressure profiles of the gas and magnetic components are given in
figure 7. The (total) gas pressure is shown in solid lines, and the magnetic
pressure - defined as $P_{\mathrm{B}}=B^{2}/8\pi$ \- is shown with dashed
lines. The red lines refer to the HMG-B simulation and the blue lines to the
HMG-snB simulation. Initially, the gas pressure dominates over the magnetic
pressure since the timescale of magnetic field amplification on this system is
large (of the order of Gyr). In both models (HMG-B and HMG-snB), the gas
pressure has a similar behaviour with higher values in the centre and
declining profile towards the outskirts of the galactic disk. The magnetic
field pressure is relatively low at the beginning of the simulation (see the
left panel in figure 7) due to the low values of the magnetic field strength,
but it rises in the outer parts of the galaxy because the magnetic field is
amplified there faster due to rotation. In the case of HMG-snB, the magnetic
pressure is lower, compared to the one in HMG-B, in the beginning since we
start with no magnetic field and seed it gradually through the stellar
feedback.

At 2 and 3 Gyr, there is a pronounced peak in the magnetic and gas pressure.
During the evolution of the galaxy, strong magnetic fields are built up in the
centre of the galaxy. This results in a magnetically driven outflow that has
been extensively discussed in Steinwandel et al. (2019) and Steinwandel et al.
(2020). This peak in the pressure components explains the peak in the star
formation profiles as our model predicts (see equation 4).

###  4.3 Evolution of the star formation rate

The total star formation rate over time is shown in figure 8. The star
formation peaks (starburst) at the beginning of the simulation because of the
high cold fraction at the initial conditions (same as SH03). A fraction of the
gas is kept in the cold star-forming phase, but there is no stellar feedback
yet in the form of SNe to heat the cold gas. This results in the initial
starburst. This can be prevented by steering ISM-turbulence at the beginning
of the simulation. However, this is deliberately not done as we want to
compare as closely as possible to the results of SH03. The star formation rate
gradually decreases with time since the star formation consumes a large part
of the gas. After almost stabilising the star formation rate, at 2.1 Gyr, the
HMG-B simulation indicates a sudden drop. This can be explained by the rising
magnetic pressure in the centre, which subsequently drives an outflow that is
correlated to adiabatic compression of the magnetic field due to a bar-
formation process. This increases the magnetic pressure in the centre to a
similar order to the magnitude of the thermal pressure (see middle and right
panel of figure 7; radial pressure profiles), leading to an outflow driven by
the magnetic pressure (see Steinwandel et al. (2020) for details).

We note that we find a difference for the model HMG-SnB compared to the
previously obtained results from Steinwandel et al. (2019). While Steinwandel
et al. (2019) reported that there is no difference in the magnetically driven
outflow between the models HMG-B and HMG-snB, we find that the outflow is
delayed by 300 Myr in the HMG-snB model compared to the HMG-B model. This is
driven by the modified star formation recipe. In the model HMG-snB, there is
no primordial magnetic field; thus, the field amplification requires longer to
build a strong enough field to drive the outflow. Additionally, the
modification in the star formation model changes the timing of the episodes of
star formation and leads to an overall slightly different evolution.

After $\sim$3 Gyr or evolution, the star formation rate for both models
reaches an asymptotic value of $\sim$1 $\mathrm{M_{\odot}/yr}$, which is in
line, within uncertainties, with observations of the star formation rate of
the Milky Way (e.g., Robitaille & Whitney, 2010; Elia et al., 2022; Zari et
al., 2023).

###  4.4 Comparison with SH03

Figure 9: Star formation rate over time for the model HMG-snB, the same
simulation run without magnetic fields (HMG-noB) and the simulations mentioned
above using the SH03 model for star formation and feedback.

In figure 9, we show the star formation rate over time for the model HMG-snB,
the same simulation run without magnetic fields (HMG-noB) and their
counterparts using the SH03 model. All simulations use the same initial
conditions described in section 3. As already mentioned, the peak of star
formation rate at the beginning of each simulation is due to the fact that
initially, a lot of gas is initiated in the cold phase, which is star-forming.
During the first Gyr of evolution, the star formation rate remains lower for
our model compared to the star formation rate for the SH03 model. This can be
explained by the fact that our model includes the molecular phase of the ISM,
which gives rise to star formation. In the simulations without magnetic fields
(orange lines, i.e., HMG-noB and HMG-noB-SH), the star formation rate stays
consistently lower due to the reasons mentioned above. The simulations,
including magnetic fields (via SN seeding), reach the same asymptotic value
after 3 Gyr of evolution. We also notice that in our model (HMG-snB), there is
a continuous decrease of the star formation rate, in contrast to the sudden
drop noticed at $\sim$2.4 Gyr in the model with the SH03 star formation
recipe, which is an impact of the magnetically driven outflow (Steinwandel et
al., 2020).

###  4.5 Magnetic Fields

Figure 10: Slice of the xy galaxy plane of the HMG model with two initial
magnetic field configurations. The top panels show the simulation with the
primordial initial magnetic field, i.e., model HMG-B, at 1 Gyr (left) and 2
Gyr (right). The bottom panels show the simulation with the SN seeding
magnetic field, i.e., model HMG-snB, at 1 Gyr (left) and 2 Gyr (right). The
colour code indicates the magnetic field strength on a logarithmic scale.
Figure 11: Total magnetic field strength of the galactic field as a function
of time. The dark blue corresponds to the simulation with the primordial
magnetic field, while the light blue corresponds to the one with the SN
seeding magnetic field. Figure 12: Edge on cross-section slices of the
simulation HMG-B that show the start and evolution of the magnetically-driven
outflow. The outflow starts at 2.14 Gyr, and we show how it is manifested
until 2.49 Gyr. The colour bar (same scale in every plot) shows the magnetic
field strength in Gauss.

In the discussed simulations, we compare two different initialisation
configurations for the magnetic fields on the galactic system. The model
HMG-B, which is initialised with the primordial magnetic field of
$10^{-9}\,\mathrm{G}$ and the model HMG-snB, which is based on magnetic
seeding from SN explosions (see the introductory part in section 4).

Figure 10 shows slices of the galactic plane for the HMG-B (top panels) and
the HMG-snB (bottom panels) for two different times of evolution, i.e. 1 Gyr
(left) and 2 Gyr (right). The colour code shows the magnetic field strength
for each particle of the computational domain. In the case of the HMG-B
simulation, the magnetic field is strong enough from the beginning, while the
HMG-snB needs more time to build up a field of similar strength. Due to the
rotation of the galactic disk, the magnetic field is amplified through the
already known $\Omega$ mechanism as a part of the $\alpha\Omega$ dynamo
(Ruzmaikin et al., 1988). In the bottom panels, it is evident that the
timescale for the amplification of the $\mathbf{B}$ in the outskirts of the
galactic plane is shorter compared to the amplification in the central parts
due to turbulent motions. The galaxy (in the textitHMG-snB model) needs
substantial time to produce stars, from which a few of them will end their
lives in SN explosions and, in turn, enrich the ISM with magnetic fields and
hot gas.

The total magnetic field strength evolution across time is shown in figure 11.
The red line corresponds to the HMG-B simulation, and the blue line
corresponds to the HMG-snB. Initially, in the HMG-B model, the primordial
magnetic field configuration is amplified in just an order of magnitude and
stays constant until it is significantly boosted at $\sim$ 2 Gyr for the HMG-B
model. The model in which the magnetic fields are initiated with the SN
seeding model (HMG-snB) shows the steep rise of the magnetic field strength
followed by an exponential growth, which is characteristic of these systems
(Pakmor & Springel, 2013; Steinwandel et al., 2019; Ntormousi et al., 2020;
Pakmor et al., 2024). The sudden magnetic field boost is also shown in the
HMG-snB simulation, although it appears later in time ($\sim$ 2.2 Gyr). In
both models, the magnetic field strength saturates of orders of $\mu G$, which
is also what is measured in observations of spiral galaxies (e.g., Chyży,
2008; Krause, 2019; Heesen et al., 2023).

In figure 12, the evolution of the magnetic field around 2.1 Gyr is shown in a
series of edge-on cross-section slices of the galaxy for the HMG-B model. As
discussed above, the magnetic field is amplified in the central part of the
galaxy, and it is correlated with a peak in the thermal gas and magnetic
pressure and thus star formation rate (last panels of figure 6 and figure 7).
Due to this steep pressure gradient, the low-density but high magnetic field
strength outflow is driven. The outflow is numerically manifested as SPMHD
particles ejected from the disk towards the gas halo, contributing to the gas
halo’s magnetic field. These particles act as carriers for the magnetic field,
allowing the outflow to evolve. The geometrical difference between our
simulations and those of Steinwandel et al. (2019) is that the outflow in our
case is not accelerated to high velocities and thus remains closer to the
galaxy (max $v_{z}\sim$ 50 km/s). The results of the simulation with the SN-
seeded magnetic fields (i.e., HMG-snB model) are similar apart from the time
delay in the outflow. Our results are also in accordance with Pakmor &
Springel (2013), who found low-density but highly magnetised rising bubbles
around the galaxy at $\sim$ 2 Gyr with similar geometry. These outflows are
also observed in different simulations, but in other cases, they are claimed
to have a different origin, for instance, in Marinacci et al. (2011), where
SN-powered bubbles drive the outflows.

###  4.6 Correlation of the magnetic field and the star formation surface
density

Figure 13: Logarithm of magnetic field strength over the surface density of
the star formation rate in the galactic disk. Red symbols show the line of
sight averaged values and blue crosses the azimuthal average. The two lines
represent theoretical scaling (Schleicher & Beck, 2013) between magnetic field
strength and star formation surface density in molecular gas (magenta line)
and neutral gas (orange line). The data are calculated after 2 Gyr of
evolution for the HMG-snB model.

We investigate the scaling of the global magnetic field with the surface
density of the star formation rate in our simulated galaxy with the full star
formation and feedback model (model HMG-snB). Theoretical studies (e.g.,
Schleicher & Beck, 2013) have shown that the magnetic field correlates with
the surface density of the star formation rate as
$\mathrm{B}\sim\Sigma^{1/3}_{\mathrm{SFR}}$.

In figure 13 the magnetic field strength over the surface density of the star
formation rate is shown in logarithmic scaling. The red dots are averages in
the line of sight across the galactic disk, and the blue crosses are azimuthal
averages. We also plot the expected theoretical scaling (Schleicher & Beck,
2013), namely, the magnetic field relation with the surface density of the
neutral gas (orange line) and the magnetic field with the surface density of
the star formation rate coming from the molecular fraction. Our results are in
very good agreement with the predictions, which are also verified by
observations (e.g., Niklas & Beck, 1997; Krause, 2015; Manna & Roy, 2023b),
which find that the total magnetic field correlates with the gas density as
$B_{\mathrm{tot}}\sim\rho^{0.5}$. It is interesting to notice that our model
tends to follow the prediction of the magenta line more closely, i.e., the
correlation between the global magnetic field and the $\Sigma_{\mathrm{SFR}}$
from molecular gas. This reflects directly on the skeleton of our model since
the star formation is driven by the molecular gas of the ISM.

###  4.7 Comparison to observations

Figure 14: Logarithm of the magnetic field strength as a function of the star
formation rate (left panel). Observational data are the galaxies from the
KINGFISH sample (Tabatabaei et al., 2017), blue dots are the model HMG-snB for
times between 1.9 and 2.2 Gyr. The right plot shows a comparison of particle
base at 2.2 Gyr, highlighting the negative correlation at a very high star-
formation rate. This is qualitatively in line with observation in NGC 1097
(Tabatabaei et al., 2018), although our simulations do not resolve the very
high star-formation density as in the observed giant molecular clouds. In
addition, the purple data points are the KINGFISH sample, divided by the
optical size, to obtain a star-formation density in the spirit of what was
done in Tabatabaei et al. (2017). The orange cross (right panel) is what we
get for model HMG-snB when we divide the total star-formation rate by the size
of the simulated galaxy.

To compute the observed mean magnetic field from the simulated galaxy (model
HMG-snB), we computed the averaged magnetic field within the stellar body of
the galaxy, i.e., taking a cylinder with a radius of $10\,\mathrm{kpc}$ and a
height of $\pm 0.5\,\mathrm{kpc}$. The mean magnetic field within our
simulated galaxy thereby compares very well with observations, as seen in the
left panel of figure 14, where we plot the mean magnetic field as a function
of global star-formation rate. Complementary to section 4.6 and figure 13, it
is shown here how the global magnetic field of the galaxy correlates with the
surface density of star-forming gas.

To confirm the observational indication that locally, in highly star-forming
molecular clouds, the local star-formation density is anti-correlated with the
local magnetic fields, we computed an equivalent star-formation surface
density for each of the star-forming SPMHD particles by dividing the star-
formation rate of each particle by the 2D area corresponding to its volume. A
comparison of this equivalent star-formation surface density with the local
magnetic field carried by the SPMHD particle is shown in the right panel of
figure 14, where we also added the measurements from NGC 1097 (Tabatabaei et
al., 2018). It is clear to see that the simulation also shows a clear anti-
correlation between the local magnetic field and the local star formation rate
at high star-formation rates. Even though the resolution of the simulation
neither allows to resolve of individual giant molecular clouds nor the
simulated galaxy to have such strong star formation as the nuclear ring in NGC
1097, it is visible that magnetic fields significantly influence and regulate
the star formation. In the spirit of Tabatabaei et al. (2017), where they
computed a global star-formation rate density by dividing through the optical
size of the galaxy, we added the KINGFISH sample (purple data-points in our
figure 14, see Figure 8 in Tabatabaei et al. (2017)) as well as the value
obtained for our simulated galaxy (orange cross in right panel of figure 14),
again showing very nice agreement.

##  5 Summary and Conclusions

We present an updated model for star formation and feedback in the
cosmological SPMHD code gadget that accounts for the magnetic field effects on
the star formation process in galaxies. We modify the model equations of the
widely used sub-grid model for star formation and feedback by Springel &
Hernquist (2003) to include non-thermal components of the ISM in the modelling
of unresolved scales. The novel feature of our model is that we include an
extra component of the multi-phase ISM apart from the hot and the cold phase,
i.e., the molecular component from which the star formation arises. Using an
observationally derived relation between the fraction of the molecular gas in
a galaxy and the total pressure in the galaxy, we include the effects of
magnetic fields via the magnetic pressure of the galaxy as a regulation of the
star formation process. This pressure-oriented approach allows the inclusion
of non-thermal components of the ISM in the star formation process. In
addition, our model is coupled to the magnetic SN seeding model (see equation
35), which deposits the magnetic field in a galaxy via SN explosions and does
not assume the primordial origin of the magnetic field in galaxies. The
effective number of SN explosions, that determine the magnetic field input, is
directly calculated from our model. Thus, our model incorporates a self-
consistent interplay between magnetic fields and star formation in galaxies.

We implemented the model in our developer’s version code gadget and tested it
by computing an isolated Milky Way-type galaxy embedded in a hot CGM
environment. Our results for the magnetic field growth and star formation rate
over time agree with previous studies of such simulations. Apart from that,
our model is able to predict both local and global properties of the magnetic
field configuration within the galaxy. Our main conclusions are summarized as
follows:

  * •

Our new star formation recipe, tested in an idealized galactic environment,
reproduces the basic galactic properties compared to previous studies (e.g.,
Pakmor & Springel, 2013; Steinwandel et al., 2019; Ntormousi et al., 2020;
Wibking & Krumholz, 2023b; Pakmor et al., 2024). Galactic magnetic fields are
amplified and reach observable values of $\mu$G. When the magnetic field is
sufficiently amplified in the galactic centre, it can drive highly magnetised
and low-density outflows, enriching the CGM with magnetic fields.

  * •

We find that the total magnetic field of the galaxy scales with the surface
density of the star formation rate, which is directly related to the gas
density of the galaxy. Our results show good agreement with theoretical
predictions (e.g., Schleicher & Beck, 2013) and observations (e.g., Tabatabaei
et al., 2017).

  * •

We observe that the local magnetic field strength is anti-correlated with the
star formation rate in the high star-formation rate regime. This indicates
that strong magnetic fields can quench star formation in giant molecular
clouds, as suggested by observations Tabatabaei et al. (2018).

Our pressure-based star formation and feedback model can be further used to
include different components of non-thermal pressure of the ISM in order to
consider all the significant processes affecting star formation. Future work
includes verifying the model in simulations of galaxies with a cosmological
background and in large-scale cosmological simulations.

## Acknowledgements

UPS is supported by a Flatiron Research Fellowship at the Center for
Computational Astrophysics of the Flatiron Institute. The Flatiron Institute
is supported by the Simons Foundation. KD acknowledges support by the COMPLEX
project from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation program grant agreement ERC-2019-AdG
882679 as well as support by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy - EXC-2094 -
390783311. MV is supported by the Fondazione ICSC National Recovery and
Resilience Plan (PNRR), Project ID CN-00000013 "Italian Research Center on
High-Performance Computing, Big Data and Quantum Computing" funded by MUR -
Next Generation EU. MV also acknowledges partial financial support from the
INFN Indark Grant, and by the Alexander von Humboldt Stiftung during the
initial phase of this work.

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Agertz & Kravtsov (2016) Agertz O., Kravtsov A. V., 2016, ApJ, 824, 79
  * Aumer et al. (2013) Aumer M., White S. D. M., Naab T., Scannapieco C., 2013, MNRAS, 434, 3142
  * Basu & Roy (2013) Basu A., Roy S., 2013, MNRAS, 433, 1675
  * Beck (2015) Beck R., 2015, A&ARv, 24, 4
  * Beck & Wielebinski (2013) Beck R., Wielebinski R., 2013, Magnetic Fields in Galaxies. p. 641, doi:10.1007/978-94-007-5612-0_13
  * Beck et al. (2012) Beck A. M., Lesch H., Dolag K., Kotarba H., Geng A., Stasyszyn F. A., 2012, MNRAS, 422, 2152
  * Beck et al. (2013) Beck A. M., Dolag K., Lesch H., Kronberg P. P., 2013, MNRAS, 435, 3575
  * Beck et al. (2016) Beck A. M., et al., 2016, MNRAS, 455, 2110
  * Bieri et al. (2023) Bieri R., Naab T., Geen S., Coles J. P., Pakmor R., Walch S., 2023, MNRAS, 523, 6336
  * Blitz & Rosolowsky (2006) Blitz L., Rosolowsky E., 2006, ApJ, 650, 933
  * Butsky et al. (2017) Butsky I., Zrake J., Kim J.-h., Yang H.-I., Abel T., 2017, ApJ, 843, 113
  * Cen & Ostriker (1992) Cen R., Ostriker J., 1992, ApJ, 393, 22
  * Chyży (2008) Chyży K. T., 2008, A&A, 482, 755
  * Crain & van de Voort (2023) Crain R. A., van de Voort F., 2023, ARA&A, 61, 473
  * Dolag & Stasyszyn (2009) Dolag K., Stasyszyn F., 2009, MNRAS, 398, 1678
  * Donnert (2014) Donnert J. M. F., 2014, MNRAS, 438, 1971
  * Draine (2011) Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium
  * Dubois & Teyssier (2010) Dubois Y., Teyssier R., 2010, A&A, 523, A72
  * Ebagezio et al. (2023) Ebagezio S., Seifried D., Walch S., Nürnberger P. C., Rathjen T. E., Naab T., 2023, MNRAS, 525, 5631
  * Elia et al. (2022) Elia D., et al., 2022, ApJ, 941, 162
  * Elmegreen & Scalo (2004) Elmegreen B. G., Scalo J., 2004, ARA&A, 42, 211
  * Federrath & Klessen (2012) Federrath C., Klessen R. S., 2012, ApJ, 761, 156
  * Gatto et al. (2015) Gatto A., et al., 2015, MNRAS, 449, 1057
  * Giammaria et al. (2021) Giammaria M., Spagna A., Lattanzi M. G., Murante G., Re Fiorentin P., Valentini M., 2021, MNRAS, 502, 2251
  * Girma & Teyssier (2024) Girma E., Teyssier R., 2024, MNRAS, 527, 6779
  * Governato et al. (2008) Governato F., Mayer L., Brook C., 2008, in Funes J. G., Corsini E. M., eds, Astronomical Society of the Pacific Conference Series Vol. 396, Formation and Evolution of Galaxy Disks. p. 453 (arXiv:0801.1707)
  * Grand et al. (2017) Grand R. J. J., et al., 2017, MNRAS, 467, 179
  * Guedes et al. (2011) Guedes J., Callegari S., Madau P., Mayer L., 2011, ApJ, 742, 76
  * Haid et al. (2016) Haid S., Walch S., Naab T., Seifried D., Mackey J., Gatto A., 2016, MNRAS, 460, 2962
  * Heesen et al. (2023) Heesen V., et al., 2023, A&A, 669, A8
  * Hernquist (1993) Hernquist L., 1993, ApJS, 86, 389
  * Hirashima et al. (2023) Hirashima K., Moriwaki K., Fujii M. S., Hirai Y., Saitoh T. R., Makino J., 2023, MNRAS, 526, 4054
  * Hirschmann et al. (2014) Hirschmann M., Dolag K., Saro A., Bachmann L., Borgani S., Burkert A., 2014, MNRAS, 442, 2304
  * Hopkins (2015) Hopkins P. F., 2015, MNRAS, 450, 53
  * Hopkins et al. (2012) Hopkins P. F., Quataert E., Murray N., 2012, MNRAS, 421, 3488
  * 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. (2022) Hopkins P. F., Butsky I. S., Panopoulou G. V., Ji S., Quataert E., Faucher-Giguère C.-A., Kereš D., 2022, MNRAS,
  * Hopkins et al. (2023) Hopkins P. F., et al., 2023, MNRAS, 519, 3154
  * Hu et al. (2017) Hu C.-Y., Naab T., Glover S. C. O., Walch S., Clark P. C., 2017, MNRAS, 471, 2151
  * Katz & Gunn (1991) Katz N., Gunn J. E., 1991, ApJ, 377, 365
  * Katz et al. (1996) Katz N., Weinberg D. H., Hernquist L., 1996, ApJ, 105, 19
  * Kauffmann et al. (1994) Kauffmann G., Guiderdoni B., White S. D. M., 1994, MNRAS, 267, 981
  * Kennicutt (1998) Kennicutt Robert C. J., 1998, ApJ, 498, 541
  * Kim & Ostriker (2015) Kim C.-G., Ostriker E. C., 2015, ApJ, 802, 99
  * Krause (2015) Krause M., 2015, Highlights of Astronomy, 16, 399
  * Krause (2019) Krause M., 2019, Galaxies, 7, 54
  * Kretschmer & Teyssier (2020) Kretschmer M., Teyssier R., 2020, MNRAS, 492, 1385
  * Krumholz & McKee (2005) Krumholz M. R., McKee C. F., 2005, ApJ, 630, 250
  * Lacey & Silk (1991) Lacey C., Silk J., 1991, ApJ, 381, 14
  * Lucas et al. (2020) Lucas W. E., Bonnell I. A., Dale J. E., 2020, MNRAS, 493, 4700
  * Manna & Roy (2023a) Manna S., Roy S., 2023a, ApJ, 944, 86
  * Manna & Roy (2023b) Manna S., Roy S., 2023b, ApJ, 944, 86
  * Marinacci et al. (2011) Marinacci F., Fraternali F., Nipoti C., Binney J., Ciotti L., Londrillo P., 2011, MNRAS, 415, 1534
  * Marinacci et al. (2014) Marinacci F., Pakmor R., Springel V., 2014, MNRAS, 437, 1750
  * Marinacci et al. (2018) Marinacci F., et al., 2018, MNRAS, 480, 5113
  * Martin-Alvarez et al. (2018) Martin-Alvarez S., Devriendt J., Slyz A., Teyssier R., 2018, MNRAS, 479, 3343
  * Martin-Alvarez et al. (2021) Martin-Alvarez S., Katz H., Sijacki D., Devriendt J., Slyz A., 2021, MNRAS, 504, 2517
  * Martizzi et al. (2015) Martizzi D., Faucher-Giguère C.-A., Quataert E., 2015, MNRAS, 450, 504
  * Mastropietro & Burkert (2008) Mastropietro C., Burkert A., 2008, MNRAS, 389, 967
  * McKee & Ostriker (1977) McKee C. F., Ostriker J. P., 1977, ApJ, 218, 148
  * Mihos & Hernquist (1996) Mihos J. C., Hernquist L., 1996, ApJ, 464, 641
  * Moster et al. (2010) Moster B. P., Macciò A. V., Somerville R. S., Johansson P. H., Naab T., 2010, MNRAS, 403, 1009
  * Murante et al. (2010) Murante G., Monaco P., Giovalli M., Borgani S., Diaferio A., 2010, MNRAS, 405, 1491
  * Murante et al. (2015) Murante G., Monaco P., Borgani S., Tornatore L., Dolag K., Goz D., 2015, MNRAS, 447, 178
  * Naab & Ostriker (2017) Naab T., Ostriker J. P., 2017, ARA&A, 55, 59
  * Navarro & White (1993) Navarro J. F., White S. D. M., 1993, MNRAS, 265, 271
  * Navarro & White (1994) Navarro J. F., White S. D. M., 1994, MNRAS, 267, 401
  * Niklas & Beck (1997) Niklas S., Beck R., 1997, A&A, 320, 54
  * Ntormousi et al. (2020) Ntormousi E., Tassis K., Del Sordo F., Fragkoudi F., Pakmor R., 2020, A&A, 641, A165
  * Oser et al. (2010) Oser L., Ostriker J. P., Naab T., Johansson P. H., Burkert A., 2010, ApJ, 725, 2312
  * Pakmor & Springel (2013) Pakmor R., Springel V., 2013, MNRAS, 432, 176
  * Pakmor et al. (2017) Pakmor R., et al., 2017, MNRAS, 469, 3185
  * Pakmor et al. (2024) Pakmor R., et al., 2024, MNRAS, 528, 2308
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 473, 4077
  * Ramesh et al. (2023) Ramesh R., Nelson D., Heesen V., Brüggen M., 2023, MNRAS, 526, 5483
  * Rieder & Teyssier (2016) Rieder M., Teyssier R., 2016, MNRAS, 457, 1722
  * Rieder & Teyssier (2017a) Rieder M., Teyssier R., 2017a, MNRAS, 471, 2674
  * Rieder & Teyssier (2017b) Rieder M., Teyssier R., 2017b, MNRAS, 472, 4368
  * Robitaille & Whitney (2010) Robitaille T. P., Whitney B. A., 2010, ApJ, 710, L11
  * Ruzmaikin et al. (1988) Ruzmaikin A., Sokolov D., Shukurov A., 1988, Nature, 336, 341
  * Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  * Scalo & Elmegreen (2004) Scalo J., Elmegreen B. G., 2004, ARA&A, 42, 275
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Schleicher & Beck (2013) Schleicher D. R. G., Beck R., 2013, A&A, 556, A142
  * Schmidt (1959) Schmidt M., 1959, ApJ, 129, 243
  * Semenov et al. (2018) Semenov V. A., Kravtsov A. V., Gnedin N. Y., 2018, ApJ, 861, 4
  * Somerville & Davé (2015) Somerville R. S., Davé R., 2015, ARA&A, 53, 51
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Springel & Hernquist (2003) Springel V., Hernquist L., 2003, MNRAS, 339, 289
  * Springel et al. (2001) Springel V., Yoshida N., White S. D. M., 2001, New Astron., 6, 79
  * Steinmetz & Muller (1995) Steinmetz M., Muller E., 1995, MNRAS, 276, 549
  * Steinwandel et al. (2019) Steinwandel U. P., Beck M. C., Arth A., Dolag K., Moster B. P., Nielaba P., 2019, MNRAS, 483, 1008
  * Steinwandel et al. (2020) Steinwandel U. P., Dolag K., Lesch H., Moster B. P., Burkert A., Prieto A., 2020, MNRAS,
  * Su et al. (2017) Su K.-Y., Hopkins P. F., Hayward C. C., Faucher-Giguère C.-A., Kereš D., Ma X., Robles V. H., 2017, MNRAS, 471, 144
  * Tabatabaei et al. (2008) Tabatabaei F. S., Krause M., Fletcher A., Beck R., 2008, A&A, 490, 1005
  * Tabatabaei et al. (2017) Tabatabaei F. S., et al., 2017, ApJ, 836, 185
  * Tabatabaei et al. (2018) Tabatabaei F. S., Minguez P., Prieto M. A., Fernández-Ontiveros J. A., 2018, Nature Astronomy, 2, 83
  * Valentini et al. (2017) Valentini M., Murante G., Borgani S., Monaco P., Bressan A., Beck A. M., 2017, MNRAS, 470, 3167
  * Valentini et al. (2018) Valentini M., Bressan A., Borgani S., Murante G., Girardi L., Tornatore L., 2018, MNRAS, 480, 722
  * Valentini et al. (2023) Valentini M., et al., 2023, MNRAS, 518, 1128
  * Vogelsberger et al. (2014) Vogelsberger M., et al., 2014, Nature, 509, 177
  * Walch & Naab (2015) Walch S., Naab T., 2015, MNRAS, 451, 2757
  * Wang & Abel (2009) Wang P., Abel T., 2009, ApJ, 696, 96
  * Wang et al. (2015) Wang L., Dutton A. A., Stinson G. S., Macciò A. V., Penzo C., Kang X., Keller B. W., Wadsley J., 2015, MNRAS, 454, 83
  * White & Rees (1978) White S. D. M., Rees M. J., 1978, MNRAS, 183, 341
  * Wibking & Krumholz (2023a) Wibking B. D., Krumholz M. R., 2023a, MNRAS, 521, 5972
  * Wibking & Krumholz (2023b) Wibking B. D., Krumholz M. R., 2023b, MNRAS, 521, 5972
  * Zari et al. (2023) Zari E., Frankel N., Rix H.-W., 2023, A&A, 669, A10