# The impact of cosmic rays on the interstellar medium and galactic outflows
of Milky Way analogues

Francisco Rodríguez Montero,1 Sergio Martin-Alvarez,2 Adrianne Slyz,1 Julien
Devriendt,1,3 Yohan Dubois,4 and Debora Sijacki5
1Sub-department of Astrophysics, University of Oxford, Keble Road, Oxford OX1
3RH, UK
2Kavli Institute for Particle Astrophysics and Cosmology (KIPAC), Stanford
University, Stanford, CA 94305, USA
3CNRS, Centre de Recherche Astrophysique de Lyon, Université de Lyon,
Université Lyon 1, ENS de Lyon, UMR 5574, F-69230 Saint-Genis-Laval, France
4Institut d’Astrophysique de Paris, UMR 7095, Sorbonne Université, CNRS, 98
bis boulevard Arago, 75014 Paris, France
5Institute of Astronomy and Kavli Institute for Cosmology, University of
Cambridge, Madingley Road, Cambridge CB3 0HA, UK
E-mail<EMAIL_ADDRESS>
(MNRAS, submitted)

###### Abstract

During the last decade, cosmological simulations have managed to reproduce
realistic and morphologically diverse galaxies, spanning the Hubble sequence.
Central to this success was a phenomenological calibration of the few included
feedback processes, whilst glossing over higher complexity baryonic physics.
This approach diminishes the predictive power of such simulations, preventing
to further our understanding of galaxy formation. To tackle this fundamental
issue, we investigate the impact of cosmic rays (CRs) and magnetic fields on
the interstellar medium (ISM) and the launching of outflows in a cosmological
zoom-in simulation of a Milky Way-like galaxy. We find that including CRs
decreases the stellar mass of the galaxy by a factor of 10 at high redshift
and $\sim 4$ at cosmic noon, leading to a stellar mass to halo mass ratio in
good agreement with abundance matching models. Such decrease is caused by two
effects: i) a reduction of cold, high-density, star-forming gas, and ii) a
larger fraction of SN events exploding at lower densities, where they have a
higher impact. SN-injected CRs produce enhanced, multi-phase galactic
outflows, which are accelerated by CR pressure gradients in the circumgalactic
medium of the galaxy. While the mass budget of these outflows is dominated by
the warm ionised gas, warm neutral and cold gas phases contribute
significantly at high redshifts. Importantly, our work shows that future JWST
observations of galaxies and their multi-phase outflows across cosmic time
have the ability to constrain the role of CRs in regulating star formation.

###### keywords:

stars: formation - galaxies: evolution - galaxies: formation - methods:
numerical - cosmic rays - galaxies: spiral - ISM: jets and outflows

††pubyear: 2023††pagerange: The impact of cosmic rays on the interstellar
medium and galactic outflows of Milky Way analogues–The impact of cosmic rays
on the interstellar medium and galactic outflows of Milky Way analogues

##  1 Introduction

Realistically modelling galaxy formation within the $\Lambda$CDM cosmological
paradigm has proven to be an extremely complex task, as its multi-scale and
multi-physics nature quickly develops a non-linear behaviour that is difficult
to grasp and model. The most sophisticated amongst modern cosmological
simulations have confronted this challenge by explicitly solving the equations
governing the evolution of matter across the vast range of scales spanning
from the cosmic large-scale structure down to the interstellar medium (ISM) of
individual galaxies. For processes such as star formation which are
unresolved, yet indispensable for galaxy formation, numerical prescriptions
are implemented. Much of their success (for a review see e.g. Vogelsberger et
al., 2020, and references therein) resides in their recognition that, without
the injection of energy by clustered supernova (SN) events, the gravitational
collapse of cold-dense gas leads to the consumption of the gas reservoir over
a free-fall time, resulting in too-efficient star formation (Kay et al., 2002;
Bournaud et al., 2010; Dobbs et al., 2011; Hopkins et al., 2011; Tasker, 2011;
Silk, 2010; Moster et al., 2013) in stark difference with observations (e.g.
Zuckerman & Evans, 1974; Williams et al., 1997; Kennicutt, 1998; Evans II,
1999; Krumholz & Tan, 2007; Evans et al., 2009). Besides the need for
reproducing the inefficient formation of stars (e.g. Stinson et al., 2006;
Teyssier et al., 2013; Springel & Hernquist, 2003), including such a
considerable energy injection through SN or active galactic nuclei (AGN)
establishes a feedback loop which is of utter importance for reproducing many
global observables of the galaxy population. Illustrative examples are the
mass-metallicity relation (Tremonti et al., 2004; Kereš et al., 2009; Mannucci
et al., 2010; Maiolino & Mannucci, 2019; Tortora et al., 2022), the metal
enrichment of the circum- (CGM) and intergalactic medium (IGM) (e.g. Veilleux
et al., 2005; Brooks et al., 2007; Chisholm et al., 2018), the quenching and
morphological transformation of late-type galaxies into early-type galaxies
(e.g. De Lucia & Blaizot, 2007; Naab et al., 2007; Dubois et al., 2013; Dubois
et al., 2016), the flattening of inner dark matter profiles (e.g. Mashchenko
et al., 2006; Governato et al., 2010; de Blok, 2010; Teyssier et al., 2013)
and the formation of bulge-less galaxies by expelling low angular momentum gas
(Governato et al., 2012).

However, this success is built on a phenomenological strategy, in which
subgrid models for feedback processes are calibrated against local observed
scaling relations (e.g. Crain et al., 2015; Pillepich et al., 2018; Rosdahl et
al., 2018; Davé et al., 2019) thus considerably curtailing predictive power.
For stellar feedback, state-of-the-art simulations are now moving into a
regime where the multi-phase structure of the interstellar medium is beginning
to be spatially resolved. This has driven the development of more physically
motivated treatments of SNe feedback (e.g. Hopkins et al., 2014; Kimm & Cen,
2014; Kimm et al., 2015), as well as other forms of stellar feedback that
could not be accounted for in the previous generations of numerical
simulations such as e.g. stellar radiation, including photoionisation heating
and radiation pressure (Murray et al., 2005; Hopkins et al., 2011; Agertz et
al., 2013; Rosdahl & Teyssier, 2015; Emerick et al., 2018; Smith et al., 2018,
2019), stellar winds (Hopkins et al., 2011; Agertz et al., 2013; Geen et al.,
2015; Fierlinger et al., 2016; Smith et al., 2019), or direct pressure from
Lyman-$\alpha$ photons (Kimm et al., 2018). Although simulations are beginning
to resolve the scales and physical processes required to model star formation
through prescriptions aimed to reproduce their formation in giant molecular
clouds, almost all aspects of the actual launching of outflows at galactic
scales and their thermodynamic properties remains a source of intense debate
in the theoretical and numerical simulations community. Furthermore, the
observations that are used to constrain stellar feedback models (e.g. Collins
& Read, 2022; Leroy et al., 2023; Thilker et al., 2023) suffer from large
limitations, with reported mass-loading factors (defined as the ratio of gas
mass ejected from the galaxy to the current rate of gas depletion by star
formation) varying by more than 2 orders of magnitude in galaxies of similar
mass (e.g. Veilleux et al., 2005; Chisholm et al., 2016). There are also
considerable difficulties involved in detecting the emission from hot ($\sim
10^{6}$ K) gas (Heckman et al., 1995; Summers et al., 2003; McQuinn et al.,
2019; Marasco et al., 2023; Concas et al., 2022), which theoretical models
argue carry a significant fraction of the kinetic energy of the wind (e.g. Kim
& Ostriker, 2018).

Cosmic rays (CRs), relativistic particles which are accelerated by diffusive
shock acceleration (e.g. Axford, 1981; Bell, 1978; Blandford & Eichler, 1987),
act as a non-thermal source of energy, and pioneering numerical investigations
(Jubelgas et al., 2008; Wadepuhl & Springel, 2011; Booth et al., 2013; Salem
et al., 2014, e.g.) of their impact indicate that they can have a significant
impact on star formation. Observations of the Milky Way and nearby galaxies
suggest that CRs and magnetic fields are in equipartition with the thermal and
turbulent energies in their mid-plane (Boulares & Cox, 1990; Basu & Roy, 2013;
Beck, 2016; Zweibel, 2017; Lopez-Rodriguez et al., 2021). Furthermore, the
peak of their spectrum at $\sim 1$ GeV leads them to exhibit smooth spatial
gradients on $\gtrsim 1$ kpc scale lengths (e.g. Ensslin, 2004). Consequently,
for some time now they have been put forward as a fundamental player in the
dynamics of extra-planar gas and the launching of outflows (e.g. Ipavich,
1975; Breitschwerdt et al., 1991; Ptuskin et al., 1997; Socrates et al., 2006;
Everett et al., 2008; Mao & Ostriker, 2018; Quataert et al., 2021b, a). These
$\sim 1$ GeV CRs behave as a relativistic fluid with a soft equation of state
(i.e. adiabatic index $\gamma_{\rm CR}=4/3$), thus making the mixture of CRs
and non-relativistic gas more compressible and undergoing smaller energetic
losses under adiabatic expansion. Additionally, for average low-redshift ISM
conditions, CRs experience longer cooling times ($\sim 10^{6-7}$ Myr) via
Coulomb and hadronic interactions (Ensslin, 2004; Schlickeiser et al., 2009)
than the radiative cooling of thermal gas. They rapidly diffuse away from
their injection sites due to diffusion lengths larger than typical molecular
clouds ($\sim 100$ pc), reaching the more dilute gas immediately above the
galactic disk, which is easier to accelerate against the gravitational
potential. Therefore, in the past few years there has been an increased effort
in modelling CR feedback in numerical studies of galactic outflows, showing
broadly consistent trends across hydrodynamic schemes and numerical
implementations of the CR Fokker-Planck equation (Skilling, 1971; Hanasz et
al., 2021, and references therein): (i) CRs thicken the gaseous disk by
altering the pressure-gravity balance, changing the morphology of the galaxy
(Booth et al., 2013; Salem et al., 2014; Ruszkowski et al., 2017; Dashyan &
Dubois, 2020; Chan et al., 2022; Farcy et al., 2022; Nuñez-Castiñeyra et al.,
2022; Buck et al., 2020; Martin-Alvarez et al., 2022). In addition, CR
feedback significantly affects the star formation history of galaxies by (ii)
reducing star formation rates (e.g. Jubelgas et al., 2008; Booth et al., 2013;
Pakmor et al., 2016; Chan et al., 2019; Dashyan & Dubois, 2020; Hopkins et
al., 2020) and by (iii) launching more heavily mass-loaded outflows, these
being denser and colder (e.g. Booth et al., 2013; Hanasz et al., 2013; Hopkins
et al., 2021a; Farcy et al., 2022; Peschken et al., 2022).

This modelling progress has established that CRs in the $\sim 1$ GeV regime
are likely an important ingredient of galaxy formation. However, a clear
understanding of its coupling with other physical processes of relevance
remains far from being explored. While isolated galaxy simulations find CRs
modelled with a constant diffusion coefficient $\kappa$ are able to affect the
star formation of galaxies over a range of masses (Booth et al., 2013;
Ruszkowski et al., 2017; Dashyan & Dubois, 2020; Farcy et al., 2022; Nuñez-
Castiñeyra et al., 2022), their inclusion in cosmological simulations shows
(e.g. Jubelgas et al., 2008; Buck et al., 2020; Martin-Alvarez et al., 2022)
that their effect on the final stellar mass of a galaxy is usually small,
preferentially affecting the most massive galaxies at low redshift (e.g.
Hopkins et al., 2020). This is likely due to complex interplay between CR-
driven outflows and cosmic gas inflows in a realistic CGM, that can be
captured only by fully self-consistent cosmological simulations (see e.g. Ji
et al., 2020, 2021; Beckmann et al., 2022; Martin-Alvarez et al., 2022).

Although $\sim 1$ GeV transport is usually modelled in the streaming
instability limit (i.e. CR-induced magnetic turbulence of the order of the
particle gyroradius), more physically motivated descriptions of CR transport
suggest that in the warm and cold phases Alfvén waves are efficiently damped
by ion and neutral collisions and by MHD turbulence (e.g. Farmer & Goldreich,
2003; Xu & Lazarian, 2022), giving rise to a variation of up to 5 orders of
magnitude in the CR mean free path (Armillotta et al., 2022) and to changes in
the acceleration of cold clouds in the CGM (Brüggen & Scannapieco, 2020;
Bustard et al., 2021). Moreover, CRs can affect the evolution of magneto-
thermal instabilities in galactic discs depending on the transport
prescription and diffusion constant (e.g. Wagner et al., 2005; Shadmehri,
2009; Kuwabara & Ko, 2020), strongly altering its structure (Commerçon et al.,
2019; Dashyan & Dubois, 2020; Nuñez-Castiñeyra et al., 2022; Martin-Alvarez et
al., 2022) and the conditions in which star formation and stellar feedback
take place. Hence, an accurate modelling of the ISM is fundamental for a
better understanding of CR feedback together with the launching and
acceleration of CR-driven multi-phase outflows (Girichidis et al., 2016; Kim &
Ostriker, 2018; Salem et al., 2014; Hopkins et al., 2021a; Martin-Alvarez et
al., 2022), which remains to be carefully studied in a cosmological context.

This is the first paper in a series in which we study the effects of CR
feedback in galaxies with a resolved multi-phase ISM using high-resolution
cosmological zoom-in simulations, with the aim to further our understanding of
the role CRs play in shaping galaxy properties. In Section 2 we outline our
Milky Way-like galaxy (called nut) formation simulation set-up, and our CR and
magnetohydrodynamics modelling. We then discuss the global effects of magnetic
fields and CRs in the evolution of nut in Section 3.1, demonstrating how the
CRs reduce the stellar mass by decreasing the star formation efficiency in
Section 3.2 and increasing the SN feedback efficiency in Section 3.3 via
modifications to the outflow phases and their more efficient acceleration. In
Section 4 we summarise our main results, discuss their limitations and the
prospect for future studies of CR feedback.

##  2 Numerical Methods

Figure 1: Overview of the CRMHD model of the nut simulation at $z=1.5$.
Background: Colour composite of the zoom-in region centred on the nut galaxy,
with gas density (silver blue), gas temperature (orange) and CR energy density
(yellow), showing 600 kpc along the x-axis. Three cold filaments are feeding
the growth of the central galaxy, where strong feedback events are driving
large-scale outflows of hot and CR-dominated gas. The effects of feedback are
also observed in smaller objects embedded in the infalling filaments. Left
inset: RGB colour composite zoom into the immediate galactic region as would
be observed with the F200W, F150W and F090W filters from JWST. This mock
observation is generated using the SKIRT radiation transfer code by assuming a
fixed dust-to-metal ratio of 0.4 and the BARE-GR-S dust model from Zubko et
al. (2004). Right inset: Edge-on view of the nut galaxy presenting different
gas quantities modelled in the CRMHD simulation. In clockwise order from top-
left corner: gas density, gas temperature, CR energy density, and magnetic
field energy density. The cold and dense disk is shown in this projection,
combined with a large hot outflow along the axis of rotation of the galaxy. CR
and magnetic energies fill the halo and are shaped by a recent outflows.

###  2.1 The RAMSES code

The full suite of new simulations presented in this work has been generated
using the ramses code111https://bitbucket.org/rteyssie/ramses/src/master/
(Teyssier, 2002), a grid-based Eulerian code with an octree AMR grid. ramses
uses collisionless particles to model the evolution of the dark matter and
stellar components. These are coupled to each other and to the gas component
through gravity described by the Poisson equation solved on the adaptive mesh.
In order to evolve in time the magneto-hydrodynamical quantities of the gas,
it makes use of a second-order Godunov scheme with the Harten-Lax-van Lear-
Discontinuity Riemann solver from Miyoshi & Kusano (2005), and the MinMod
slope limiter to reconstruct the cell-centred properties. The implementation
of magnetic fields in ramses uses the constrained transport (CT) method, which
models magnetic field values as cell face-centred quantities (Teyssier et al.,
2006; Fromang et al., 2006). This method ensures that the solenoidal condition
(i.e. $\vec{\nabla}\cdot\vec{B}=0$) is fulfilled to within numerical precision
errors, preventing the numerically spurious growth of magnetic field and cell
noise (see Tóth et al., 2000; Mocz et al., 2016, for a comparison of CT with
other MHD methods).

####  2.1.1 Cosmic ray magneto-hydrodynamics

We make use of the implementation of CRs from Dubois & Commerçon (2016) and
Dubois et al. (2019), in which CRs are treated as a separate fluid represented
by one additional energy equation in the set of MHD fluid equations:

$\displaystyle\frac{\partial\rho}{\partial
t}+\vec{\nabla}\cdot(\rho\vec{u})=0,$ (1)
$\displaystyle\frac{\partial\rho\vec{u}}{\partial
t}+\vec{\nabla}\cdot\left(\rho\vec{u}\vec{u}+P_{\rm
tot}-\frac{\vec{B}\vec{B}}{4\pi}\right)=\rho\vec{g},$ (2)
$\displaystyle\frac{\partial e}{\partial
t}+\vec{\nabla}\cdot\left((e+P)\vec{u}-\frac{\vec{B}(\vec{B}\cdot\vec{u})}{4\pi}\right)$
$\displaystyle=\rho\vec{u}\cdot\vec{g}-P_{\text{CR}}\vec{\nabla}\cdot\vec{u}-\vec{\nabla}\cdot\vec{F}_{\text{CR,d}}+\mathcal{L}_{\text{rad}},$
(3) $\displaystyle\frac{\partial\vec{B}}{\partial
t}-\vec{\nabla}\times(\vec{u}\times\vec{B})=0,$ (4)
$\displaystyle\frac{\partial e_{\text{CR}}}{\partial
t}+\vec{\nabla}\cdot(e_{\text{CR}}\vec{u}+(e_{\text{CR}}+P_{\text{CR}})\vec{u}_{\text{st}})$
$\displaystyle=-P_{\text{CR}}\vec{\nabla}\cdot\vec{u}-\vec{\nabla}\cdot\vec{F}_{\text{CR,d}}+\mathcal{L}_{\text{st}}+\mathcal{L}_{\text{rad,CR},}$
(5)

where $\rho$ is the gas mass density, $\vec{u}$ is the gas velocity vector,
$\vec{g}$ the gravitational field,
$\vec{u}_{\text{st}}=-\vec{u}_{\text{A}}\text{sign}(\vec{b}\cdot\vec{\nabla}e_{\rm
CR})$ is the streaming velocity (where $\vec{u}_{\text{A}}$ is the Alfvén
velocity and $\vec{b}=\vec{B}/|\vec{B}|$ the magnetic unit vector), $\vec{B}$
is the magnetic field,
$e=0.5\rho|\vec{u}|^{2}+e_{\text{th}}+e_{\text{CR}}+|\vec{B}|^{2}/(8\pi)$ is
the total energy density, $e_{\text{th}}$ is the thermal energy density, and
$e_{\text{CR}}$ is the CR energy density, $P_{\rm
tot}=P_{\text{th}}+P_{\text{CR}}+P_{\text{mag}}$ is the sum of the thermal
$P_{\text{th}}=(\gamma-1)e_{\text{th}}$, CR
$P_{\text{CR}}=(\gamma_{\text{CR}}-1)e_{\text{CR}}$ and magnetic
$P_{\text{mag}}=|\vec{B}|^{2}/(8\pi)$ pressures, with adiabatic indexes
$\gamma=5/3$ and $\gamma_{\rm CR}=4/3$ for an ideal non-relativistic and
relativistic gas, respectively. This modelling approximation can be adopted
for the evolution of low-energy ($\sim$GeV) CRs, as their interaction with the
thermal gas is mostly collisionless and mediated by pitch-angle scattering by
the local magnetic field. A practical way of interpreting the terms on the
right side of Equations (3) and (5) is to treat them as source terms, with
$P_{\text{CR}}\vec{\nabla}\cdot\vec{u}$ the CR pressure work, $\vec{F}_{\rm
CR,d}=-\kappa\vec{b}(\vec{b}\cdot\vec{\nabla}e_{\rm CR})$ the anisotropic
diffusion flux ($\kappa$ is the diffusion coefficient, assumed constant in
this work), $\mathcal{L}_{\rm rad}=\mathcal{L}_{\rm rad,th}+\mathcal{L}_{\rm
rad,CR\rightarrow th}$ the total radiative losses which includes thermal and
CR radiative losses, respectively222The modelling of radiative losses is
described in Section 2.2.. The CR loss term encompasses both energy losses due
to hadronic and Coulomb interactions $\mathcal{L}_{\rm CR}=-7.51\times
10^{-16}(n_{\rm e}/\text{cm}^{-3})(e_{\rm CR}/\text{erg cm}^{-3})\text{ erg
s}^{-1}$ (see e.g., Ensslin, 2004; Guo & Oh, 2008) transferred to the thermal
component at a rate $\mathcal{H}_{\rm CR\rightarrow th}=2.63\times
10^{-16}(n_{\rm e}/\text{cm}^{-3})(e_{\rm CR}/\text{erg cm}^{-3})\text{ erg
s}^{-1}$.

In the CR energy equation (Equation (5)), $\vec{F}_{\rm st}=(e_{\rm CR}+P_{\rm
CR})\,\vec{u}_{\rm st}$ is the streaming advection term, which emulates the
self-confinement of CR transport by the self-excited Alfvén waves (Kulsrud &
Pearce, 1969; Skilling, 1975). As we have discussed in Section 1, several
mechanisms may play an important role in damping these self-excited Alfvén
waves in different gas phases (e.g. Armillotta et al., 2022), allowing for
streaming velocities exceeding $\vec{u}_{\rm A}$. As a preliminary step
towards a deeper understanding of CR feedback in a cosmological context, we
assume that CR streaming velocities cannot reach super-Alfvénic velocities,
and leave the implementation of self-consistent CR diffusion coupling (like
the one presented in Farber et al., 2018; Xu & Lazarian, 2022; Armillotta et
al., 2022) for future work. Furthermore, since CRs scatter on the Alfvén
waves, they experience a drag force, transferring energy to the thermal gas at
a rate $\mathcal{L}_{\rm st}=-\text{sign}(\vec{b}\cdot\vec{\nabla}e_{\rm
CR})\vec{u}_{\rm A}\cdot\vec{\nabla}P_{\rm CR}$. The typical value of
$\kappa\simeq 3\times 10^{28}(R_{\rm CR}/1\text{GV})^{0.34}\text{cm}^{2}\text{
s}^{-1}$ depends on the particle rigidity $R_{\rm CR}$ (Jóhannesson et al.,
2019), which for 1 GeV proton has a value of $R_{\rm CR}=1\text{ GV}$.
Previous efforts to constrain $\kappa$ in CR simulations by comparing Fermi
LAT $\gamma$-ray observations(e.g. Fire-2, arepo and ramses, Chan et al.,
2019; Werhahn et al., 2021; Nuñez-Castiñeyra et al., 2022, respectively) have
resulted in large discrepancies in the inferred diffusion coefficient favoured
by observations. Nuñez-Castiñeyra et al. (2022) argue that, although
differences in the assumed CR spectrum or $\gamma$-ray emissivity per atom can
reduce the discrepancy by factors of up to 2.5, these are not sufficient to
bring different simulations to agreement. They suggest that differences in the
thermal state of the gas and the magnetic field properties, achieved by
different MHD codes, could be the cause of these differences. The exploration
of how different CR transport mechanisms and diffusion coefficients affect the
$\gamma$-ray luminosity is the scope of an upcoming paper (Rodríguez Montero
et al., in prep). We adopt here a value of $3\times
10^{28}\text{cm}^{2}\,\text{s}^{-1}$ favoured by the results of isolated CR
simulations using similar versions of ramses to the one used in this work
(Dashyan & Dubois, 2020; Farcy et al., 2022; Nuñez-Castiñeyra et al., 2022),
and leave the analysis of a higher $\kappa=3\times
10^{29}\text{cm}^{2}\,\text{s}^{-1}$ (favoured by the Fire-2 simulations Chan
et al., 2019; Hopkins et al., 2021b) to future work.

###  2.2 The nut suite of simulations

The new suite of simulations (see Table 1 for details) presented in this work
uses the initial conditions of the nut suite (Powell et al. (2011) and Geen et
al. (2013)) generated at $z=499$ with cosmological parameters consistent with
the results of WMAP5 (matter density $\Omega_{\rm m}=0.258$, dark energy
density $\Omega_{\Lambda}=0.742$, baryon density $\Omega_{\rm b}=0.045$,
amplitude of the linear power spectrum at the $8\,\rm Mpc$ scale
$\sigma_{8}=0.8$ and Hubble constant $H_{0}=72\,\rm km\,s^{-1}\,Mpc^{-1}$
Dunkley et al., 2009). The computational box has a side of comoving length 9
Mpc $h^{-1}$ with a spherical region of approximately 2.7 comoving Mpc
$h^{-1}$ in diameter defining the volume in which high resolution is
performed. Within this zoom-in region a $M_{\rm vir}(z=0)\simeq 5\times
10^{11}M_{\odot}$ DM halo hosts a Milky Way-like galaxy (from this point
simply identified as the nut galaxy) with a DM particle resolution of $m_{\rm
DM}=5.5\times 10^{4}M_{\odot}$. To probe the multi-phase structure of the ISM,
we allow refinement within this region to be triggered via a quasi-Lagrangian
method, for which cells with total mass (i.e. gas, stars and dark matter)
above $8m_{\rm DM}\Omega_{\rm m}/\Omega_{\rm b}$ are subdivided into 8 equal
cells, reaching a maximum physical resolution of cell size $\Delta
x_{\max}\simeq 23$ pc approximately constant in proper units (an extra level
of refinement is triggered for expansion factors of $a_{\rm
exp}=0.1,0.2,0.4,0.8$).

In Fig. 1 we show an overview of our fiducial simulation with CRs and MHD
centred on nut. The main galaxy is fed by three large filaments (gas density
in silver), and at their intersection there is a gas halo filled with high
temperature gas (in orange) and CR energy (yellow) from recent feedback
events. Additionally, we show zoom-in insets of the galactic region of nut
seen face-on (left circle) as it would be seen through JWST F200W (red), F150W
(green) and F090W (blue) filters, and edge-on (right circle) density-weighted
projections of the raw simulation data. The mock JWST RGB image has been
obtained with the radiation transfer code skirt (Camps & Baes, 2015), assuming
a fixed dust-to-metal ratio of 0.4 and the BARE-GR-S dust model from Zubko et
al. (2004). At $z=1.5$ nut shows signatures of a large feedback event: hot
bipolar outflows emerge perpendicular to the cold, dense and magnetised disk.
These outflows clear denser and colder gas in the CGM, creating bubbles seen
in the gas density, magnetic energy and CR energy. The face-on view through
SDSS filters shows a late-type spiral with a nuclear cluster at its centre
dominated by an old stellar population, and a disk filled with lanes and
filaments of dense dust absorption and small regions of recent star formation.

For all our runs, we include metal-dependent radiative cooling with
interpolated cooling curves from cloudy (Ferland et al., 1998), down to $\sim
10^{4}$ K. Below this temperature, gas is allowed to cool further via metal
fine structure atomic transitions (Rosen & Bregman, 1995) (to a minimum of 15
K) and via adiabatic expansion. We employ a uniform ultraviolet background
that is turned on at a redshift of $z=9$ (Haardt & Madau, 1996), with gas
denser than $n_{\rm H}=0.01$ cm-3 considered to be self-shielded from UV
radiation with an exponential dependence on $n_{\rm H}$. Furthermore, despite
the high resolution of our simulation, we lack the cooling modelling
capabilities to resolve primordial (Population III) SN explosions in minihalos
(e.g. Whalen et al., 2008; Wise et al., 2012) that would potentially raise the
initial gas metallicity to $\sim 10^{-3}\,Z_{\odot}$ and hence increase
cooling at high redshift. To avoid this issue, we initialise the simulation
with a uniform metallicity floor of $10^{-3}\,Z_{\odot}$.

The formation of stellar particles from the gas only takes place in the cells
at the highest level of refinement when the gas is found to be gravitationally
unstable. We determine this following the magneto-thermo-turbulent (MTT) star
formation model presented in Kimm et al. (2017) and in Martin-Alvarez et al.
(2020) for its MHD version. Star formation can only happen where the gas
number density is above 10 cm-3 and the gravitational pull is larger than the
combined support from turbulent, thermal, magnetic and CR pressure333We have
checked the change in efficiency due to the addition of CR pressure to the
original MTT model, and found that it commonly decreases the efficiency by
$\sim 10-30$%, although in rare occasions the support of CR pressure can be
enough to turn the efficiency to 0.. We determine the rate of gas to stellar
mass conversion following a Schmidt law (Schmidt, 1959)

$\displaystyle\dot{\rho}_{\rm star}=\epsilon_{\rm ff}\frac{\rho_{\rm
gas}}{t_{\rm ff}},$ (6)

where we allow $\epsilon_{\rm ff}$ to be a local parameter that is determined
by the gas cell conditions as predicted by the multi-free-fall model from
Hennebelle & Chabrier (2011) (generalised to the model of PN11 by Federrath &
Klessen (2012)), and refer the reader interested in further details of this
model to appendix B in Martin-Alvarez et al. (2020). Ultimately, the
conversion of gas to stars is modelled by sampling stochastically a Poisson
mass-probability distribution for the final stellar particle mass (Rasera &
Teyssier, 2006), with the minimum stellar mass at creation fixed at $4.5\times
10^{3}M_{\odot}$.

As discussed in the introduction, instead of performing a calibration process
to determine artificial ‘boosting’ factors to the canonical type II SN energy
of $10^{51}$ erg that may allow, e.g., the recovery of a stellar-to-halo-mass
(SMHM) relation in agreement with local observations (Rosdahl et al., 2018),
we model SN feedback using the ‘mechanical’ approach of Kimm & Cen (2014).
Considering the local gas conditions and the available spatial resolution, it
determines whether the Sedov-Taylor expansion phase of the SN can be captured
directly (Blondin et al., 1998; Thornton et al., 1998) or if numerical
overcooling will artificially affect its evolution. In the latter case,
momentum prescribed by the snowplow phase is injected. In order to determine
when a stellar particle should experience a SN feedback event, we follow the
multiple explosions method by Kimm et al. (2015), in which we use the real
delay of SNe given by the population synthesis code starburst99 (Leitherer et
al., 1999), with SNe taking place as early as 3 Myr and as late as 50 Myr.
These SN events inject mass, momentum and energy into the host and a maximum
number of 48 neighbouring cells, with a canonical efficiency of
$\varepsilon_{\rm SN}=10^{51}\text{ erg}/(10M_{\odot})$ and a fraction of 0.2
of the SN mass returned to the SN host cell. Additionally, a fraction 0.075 of
this total mass corresponds to metals, which are followed via a passive scalar
advected with the gas flow by the Godunov solver. We further assume a Kroupa
IMF (Kroupa, 2001) for the underlying stellar population of our star
particles.

The ideal MHD induction equation prohibits the formation of magnetic fields in
the absence of an initial non-zero $\vec{B}$ field. Several mechanisms have
been presented over the last decades as candidates to generate galactic
magnetic fields of the order of $10^{-7}-10^{-5}$ G (e.g. Mulcahy et al.,
2014; Beck, 2016). These can be classified by their cosmological (e.g.,
electroweak or quantum chromo-dynamics phase transitions Vachaspati, 1991;
Olesen, 1993) or astrophysical (plasma seeding processes like Biermann battery
Biermann, 1950; Pudritz et al., 1989) nature (see Durrer & Neronov, 2013;
Vachaspati, 2021, for recent reviews on the matter). Whether magnetic fields
in galaxies are primordial in nature or rapidly amplified by astrophysical or
dynamo processes (Pakmor & Springel, 2013; Rieder & Teyssier, 2017; Vazza et
al., 2017; Martin-Alvarez et al., 2018, 2021), the theoretical expectation is
that they reach $\sim 10^{-6}$ G strengths no later than a few hundred Myr
after galaxy formation (Martin-Alvarez et al., 2022). We permeate the
simulation box with an uniform magnetic field along the $z$ axis444In this
work we do not explore the effect of different orientations, as this does not
significantly alter the evolution of nut (Martin-Alvarez et al., 2020), but
the magnetisation of the large scale environment changes.. In order to mimic
such rapid amplification, we employ a uniform primordial magnetic field of
strength $B_{0}=3\times 10^{-12}$ G. This simple magnetic field initialisation
directly provides the expected strengths of $1-10$ $\mu$ G in the nut galaxy
shortly after its formation (Martin-Alvarez et al., 2021).

In order to properly capture CR evolution, we activate the CR solver when the
first stars are formed: shortly prior to their first injection into the
simulation. This is because SNe are the exclusive source of CRs in our
simulations. We assume that each SN event injects a fixed fraction of its
total energy in the form of CR energy into the gas cell hosting the SN
particle. We choose this fraction to be $f_{\rm CR}=0.1$, which leads to a CR
injected energy per SN event of $10^{50}$ erg. This value is in fair agreement
with the estimated acceleration efficiency of CRs by non-relativistic shocks
(e.g. Morlino & Caprioli, 2012; Caprioli & Spitkovsky, 2014). In order to keep
the energy injection for type II SN fixed to the canonical value, the energy
available for thermal energy and radial momentum injection by the mechanical
feedback prescription is reduced accordingly.

Simulation snapshots are analysed using
Ozymandias555https://github.com/Currodri/ozymandias a new cosmological
simulation analysis package which allows for a smooth interfacing of halo and
galaxy selection codes, in-depth analysis algorithms and simulated mock
observations. It is developed with an easy-access philosophy in mind,
combining Python for accessing the data catalogues and Fortran 90, openMP
parallelised routines to access the raw ramses data. It takes a cosmological
simulation snapshot and outputs a convenient galaxy and halo catalogue
following the HDF5666https://www.hdfgroup.org/solutions/hdf5/ format standard,
which makes possible the study of many physical properties in the simulation
without the need to read back again the original binary files. Additionally,
these files are usually less than 0.3% the size of the original outputs,
making simple analysis and comparison between simulations much less memory and
computationally expensive. Ozymandias has been written to work with two
different halo and galaxy finders (although the I/O for halo catalogues is
fully compatible with other halo finders, providing the data structure
information): the HaloMaker code by Tweed et al. (2009), a 3D density-
threshold structure finder extended to use the shrinking spheres centring
method (Power et al., 2003); and the VELOCIraptor code (Elahi et al., 2019;
Cañas et al., 2019) in its enhanced version by Rhee et al. (2022), a robust,
scale-free, phase-space structure-finder. Once halos and galaxies are
identified, a hierarchical structure is built between them, linking galaxies
with their host halos and linking satellites with central galaxies. Then, in
order to allow for an easy tracking of each galaxy/halo across cosmic time,
particle lists for each identified object are compared between consecutive
snapshots, allowing for an easy access to the progenitors in the previous
snapshots. We use this package to identify the main galaxy, and refer to the
galactic region as the sphere with radius $0.2R_{\rm vir,DM}$, $R_{\rm
vir,DM}$ being the host DM halo virial radius. Quantities such as the galaxy
stellar mass $M_{*}$ and gas mass $M_{\rm gas}$ are measured within this
spherical volume at every snapshot.

Table 1: Runs of our new nut suite of simulations studying CR feedback included in this work. Columns indicate from left to right the simulation name, their employed solver and the details of the CR physics if CRs are present. All simulations use the same initial conditions, a MTT star formation prescription and the canonical configuration for mechanical feedback prescription. For simulations using the MHD or the CRMHD solver, all simulations employ a uniform comoving initial magnetic field of $B_{0}=3\times 10^{-12}$ G. Simulation |  Solver |  CR physics
---|---|---
HD |  Hydro |  ✗
MHD |  MHD |  ✗
CRMHD |  CRMHD |  Streaming and $\kappa=3\times 10^{28}$ cm2 s-1

##  3 Results

###  3.1 The global effect of cosmic rays in the evolution of nut

Figure 2: Projections of the nut galaxy for the three main models explored in
this work: HD (top 2 rows), MHD (middle 2 rows), and CRMHD (bottom 2 rows),
generated at $z=1.5$. For each model, the plot includes density-weighted
projections of gas density (first column), gas temperature (second column),
and gas metallicity (third column), 20 kpc on a side. Additionally, we present
in the last column RGB synthetic Skirt observations in the JWST F200W, F150W
and F090W filters, respectively, for a zoom region of 10 kpc on a side. Face-
on and edge-on views are shown for all models. Dashed white circles indicate
the half-light radius $r_{50}$ obtained by fitting the RGB images with a 2D
Sérsic profile. For each simulation the stellar mass within $0.2\,R_{\rm
vir,DM}$ and galaxy compactness $\log{\Sigma_{1.5}}$ (Barro et al., 2013), is
listed, with a larger value indicating a more compact stellar emission.
Compared to HD and MHD simulations, the disk of nut shows much more structure
in the CRMHD simulation, dominated by a more diffuse ISM with less compact and
dense clumps. The CGM is also strongly altered, with the cold dense filaments
of HD and MHD replaced by a diffuse and cold halo more uniformly enriched by
metals. The mock observations of HD and MHD show a galaxy dominated by an old
stellar population and a large bulge, while the nut galaxy in CRMHD is
transformed into a less compact system dominated by a disk. In this case, star
forming regions are scattered across spiral arms and clumps all over the disk.

In Fig. 2 we show density-weighted projections of nut at $z=1.5$ for the three
simulations (HD, MHD and CRMHD), comparing face-on views (line-of-sight
parallel to the angular momentum vector of all baryons assigned to the galaxy
by the halo finder) and edge-on views. We have also included mock JWST RGB
images obtained with SKIRT following the same parameter configuration as for
Fig. 1. At $z=1.5$, the nut galaxy has formed a thin gas disk in all
simulations, that has large-scale spiral features with high density clumps and
low density inter-arm regions. In the HD and MHD simulations (top and middle
rows) the galaxy has the general properties of a disk galaxy within a
virialised halo: a cold ($T<10^{4}$ K) thin disk with close to solar
metallicity Z⊙ is embedded within a low density, hot ($T>10^{6}$ K) CGM. The
edge-on views aid in identifying long filaments of cold, dense and low
metallicity gas inflowing into the galaxy. In the CRMHD simulation (bottom
rows), while the dense and cold region of the disk can still be clearly seen
in the projections, there are significant differences with respect to the no-
CR simulations. Firstly, the gas disk shows much more filamentary
substructure, and in comparison to the HD and MHD simulations the maximum
value of density in the clumps is substantially reduced, as well as the volume
occupied by low density, hot gas in the inter-arm regions. Secondly, the high
density centre of the galaxy in HD and MHD extends to $\sim 1-2$ kpc, while in
CRMHD it is $3$ times smaller. However, the most striking difference to the HD
and MHD simulations is the thermodynamical state of the CGM: the CRMHD run has
no clear inflow filaments, and it is instead dominated by homogeneous and
denser gas at $T\sim 10^{5}$ K. Against this lower temperature ‘background’,
recent hot outflows can clearly be distinguished in the edge-on view, although
the lower average metallicity of the disk gas compared to HD and MHD makes the
metal enrichment of outflows less prominent. Instead, the CGM of nut is now
more uniformly enriched to $\sim 0.1$ Z⊙.

Furthermore, the mock JWST images of the HD and MHD simulations show nut as a
disk galaxy dominated by a massive bulge and an old stellar population, with
dark, dense filaments seen as dust lanes in the edge-on views. Also, massive
stellar clumps formed by large-scale disk instabilities persist due to
inefficient feedback (e.g. Inoue et al., 2016; Bournaud et al., 2013). In
comparison, the CRMHD simulation shows a dramatic bulge reduction, with the
central luminosity of the galaxy reduced by 0.5 dex. Additionally, instead of
having massive clumps within its disk, small star forming regions are seen all
across the disk, with the luminosity of young stellar population more
extended. To better quantify this important difference, we have measured the
compactness of the stellar emission of nut using the total luminosity in the
JWST filters. By fitting a 2D Sérsic profile, we obtain half-light radii
$r_{50}$ for the three face-on observations, which allow us to compute the
compactness as defined in Barro et al. (2013) and used for SDSS galaxies
(Baldry et al., 2020):

$\displaystyle\log\Sigma_{1.5}=\log(M_{*}/{\rm M}_{\odot})-1.5\log(r_{50}/{\rm
kpc})\,,$ (7)

with a larger value of $\log\Sigma_{1.5}$ corresponding to a more compact
galaxy. The HD simulation has the highest level of compactness, with the MHD
simulation being slightly less compact but very similar as their total $M_{*}$
and $r_{50}$ are comparable. When considering the CRMHD simulation, a
significantly reduced $r_{50}$ and a stellar mass three times smaller than in
the HD simulation results in the least compact of our galaxies.

Figure 3: Evolution of the stellar (left panel) and gas (right panel) masses
contained within the galactic region (i.e. $<0.2R_{\rm vir,DM}$) of the nut
galaxy in the HD (blue line), MHD (purple line) and CRMHD (green line)
simulations. The right panel distinguishes the total gas mass (solid line)
from the cold gas mass (dashed line), as defined by the constant entropy limit
of the cold phase of the ISM by Gent (2012). In all runs stellar mass grows
rapidly up to $z\sim 4$, which marks the end of the accretion dominated phase,
but the simulation with CRs has a stellar mass 1 dex lower than the HD and MHD
runs. During this rapid growth period, the total gas mass of the CRMHD run is
lower than in the other two simulations, but the largest difference is seen in
the cold gas mass, which is reduced by $\sim 1$ dex at $z\gtrsim 4$. Below
$z\sim 3$, the CRMHD run experiences a fast growth in cold gas mass, which
translates into a faster growth of stellar mass at lower redshifts, causing
the gas mass at $z=1.5$ to be similar to that in the HD and MHD runs, whilst
the stellar mass remains 4 times lower. Figure 4: Stellar mass to halo mass
function of the nut galaxy for our fiducial models: HD (blue), MHD (magenta)
and CRMHD (green) shown in 4 redshift ranges ($z\sim 8,6,3$ and 1.5 –
separated by dashed vertical lines), representative of their evolution. We
have superimposed, as a visual guide, the results from abundance matching of
Moster et al. (2013) (blue-shaded region) and Behroozi et al. (2013) (orange-
shaded region), evaluated at $z=1.5$. Coloured dashed orange and blue lines
indicate the stellar masses at which observational mass functions are subject
to incompleteness. The HD and MHD runs are able transform the majority of
their baryonic content into stars, while the CRMHD run from early on ($z\sim
8$) sits $\sim 1$ dex below, approaching the prediction by Behroozi et al.
(2013) at $z\sim 1.5$.

We begin by briefly revisiting the evolution of the HD and MHD runs, as these
have already been explored for the same initial conditions in previous works
(Powell et al., 2011; Kimm et al., 2015; Martin-Alvarez et al., 2018). The
stellar (left panel) and gas (right panel) masses inside the galactic region
are shown for the HD model with blue lines in Fig. 3. The stellar mass growth
follows the common trend for a MW-like halo (e.g. Kim & Ostriker, 2015;
Martin-Alvarez et al., 2018): 1) accretion phase: significant gas accretion at
high redshifts allows for a fast growth of stellar mass, reaching
$10^{10}\text{M}_{\odot}$ at $\sim 1.5\,\rm Gyr$ after the start of the
simulation; 2) feedback phase: feedback from SNe and the formation of a virial
shock after $z\sim 4$ slows down star formation, leading to self-regulation
with $M_{*}\sim 4\times 10^{10}\text{M}_{\odot}$ at $z=1.5$, in agreement with
the equivalent MFBm model in Kimm et al. (2015). The evolution of $M_{\rm
gas}$ is also in agreement with the findings from Kimm et al. (2015),
exhibiting a rapid growth until $z\sim 5$. Subsequently, the fast depletion
timescale reduces the gas mass inside the galactic region until $z\sim 3$. A
new inflow of gas mass takes place from $z=3$ to $z=2$, coinciding with the
last major merger that the nut galaxy experiences (Martin-Alvarez et al.,
2018, e.g.). Fig. 3 includes the MHD run with a primordial magnetic field of
$3\times 10^{-12}$ G but no cosmic rays. This illustrates that the moderate
magnetisation employed here does not significantly affect the global
properties of the galaxy (Su et al., 2017; Pillepich et al., 2018; Hopkins et
al., 2020; Martin-Alvarez et al., 2018).

However, the addition of CRs gives rise to a radical difference in the star
formation history (SFH) of our nutgalaxy: the steepness of the stellar mass
growth is reduced during the accretion phase, resulting in $M_{*}\simeq
2\times 10^{9}\text{M}_{\odot}$, $\sim 1$ order of magnitude lower with
respect to the HD and MHD models at $z\sim 4$. During this period of high-
redshift gas accretion, the CRMHD run has only a $\sim\times 2-4$ deficit in
$M_{\rm gas}$ inside the galactic region. Therefore the decrease in $M_{*}$
cannot be solely attributed to the depletion of $M_{\rm gas}$ from the
galactic region. We also include on the right panel the cold gas mass (dashed
lines) as defined by the constant entropy limit777The cold phase specific
entropy limit is given by $s\leq\qty{4.4e8}{\,{}^{-1}{}^{-1}}$, e.g.
corresponding to a density of 7.8e-23 ^-3 at 300 K. of Gent (2012), showing
that for the CRMHD simulation the cold gas mass decreases by $\sim 1$ dex
compared to HD and MHD at $z\gtrsim 4$. This points at the depletion of the
cold gas mass in the ISM as the major cause of the lower stellar mass growth
at these cosmic times. Following this period of strong star formation
suppression, the CRMHD run exhibits a similar period of gas accretion
beginning at $z\sim 3$ as the rest of the runs, which eventually translates
into a new phase of steep $M_{*}$ growth. This growth of stellar mass is
higher in relative terms for the CRMHD run than for the HD and MHD runs,
reducing the gap between stellar masses to a factor of $\sim 4$, while the
cold gas mass becomes similar in the three runs.

To provide further context for the evolution of the nut stellar masses, Fig. 4
shows the stellar mass ($M_{*}$) to halo mass ($M_{\rm halo}$) relation for
the HD (blue diamonds), MHD (purple stars) and CRMHD (green triangles) runs at
four different redshift ranges from $z=1.5$ to $z=8$. We define the halo mass
as the full mass enclosed in gas and DM within $R_{\rm vir,DM}$ and the
stellar mass of the galaxy. These redshift ranges are separated by vertical
dashed lines with their central redshift values indicated at the bottom of the
panel. Results from the abundance matching models by Behroozi et al. (2013)
(orange shaded region) and Moster et al. (2013) (blue shaded region) at
$z=1.5$ are included for reference. The HD and MHD runs show very similar
behaviour across the redshift range presented, with only a small stellar mass
suppression by magnetic fields. Both feature large baryonic conversion
efficiencies, close to the line of total conversion of gas into stars for a
given halo, given by $f_{\rm b}M_{\rm halo}$, where $f_{\rm b}$ is the
universal baryon fraction defined as $\Omega_{\rm b}/\Omega_{\rm m}$, the
ratio of the cosmological baryon and matter density parameters. This result is
consistent with the well established understanding of non-calibrated,
canonical SN feedback: it fails to reproduce the results of abundance matching
models, leading to a too massive stellar mass for a given dark matter halo
mass. The HD run begins with $\sim 32$% baryon conversion efficiency and grows
to $\sim 56$% by $z=2$, while the MHD run has a slightly lower efficiency of
$\sim 47$%. The CRMHD run shows a consistent decrease of the baryonic
conversion efficiency across all redshifts, with a maximum value of $\sim 14$%
at $z=2$ and a minimum of 5% at $z=6$. This significant decrease in stellar
mass for a given DM halo in the CRMHD model results in a much better match to
the semi-empirical constraints of Behroozi et al. (2013) by $z=1.5$. Despite
the decrease of almost an order of magnitude in $M_{*}$ for the CRMHD run
compared to the HD and MHD runs, the CRMHD model does not follow the ‘natural’
extrapolation of the abundance matching models to lower halo masses. It is
important to point out that reliable constraints on the number density of
galaxies hosted by low mass halos of $10^{10}\text{M}_{\odot}$ between $z=6$
and $z=8$ are scarce, adding uncertainty in abundance matching models at high
redshift. Additionally, the observational stellar mass functions used in these
works are subject to observational incompleteness below
$10^{7}\text{M}_{\odot}$ (shown with the dashed orange and blue lines). We
expect observations by JWST and Rubin-LSST will be able to remedy the
situation in a near future.

###  3.2 How cosmic rays alter the gas distribution and the star formation
efficiency

Figure 5: Stacked gas temperature-density phase diagrams of the nut galaxy
($<0.2R_{\rm vir,DM}$) in HD (left panel), MHD (middle panel) and CRMHD (right
panel) at $z\sim 2.5$. Colour maps indicate the total gas mass in each bin,
with lower masses shown in dark blue and higher masses in bright aquamarine.
We include a dashed white line corresponding to the median of the 2D profile.
Gas phases are separated by lines of constant specific entropy given by Gent
(2012), into hot (top left, red), warm (central, orange), and cold (bottom
right, blue). Black contours indicate the regions where $\sim 80$% of the gas
with $\epsilon_{\rm SF}^{\rm MTT}\geq 1$% lie. These show how the CRMHD run
has much lower total gas mass in regime favourable to star formation.

As shown in Section 3.1, the CRMHD simulation contains a factor of $\sim 2-3$
lower amount of gas in the galactic region of nut (compared to HD and MHD)
between $z=8$ and $z=4$. This epoch also is where we find the largest
discrepancy in $M_{*}$ between the CRMHD and HD/MHD runs. However, the 10
times smaller stellar mass measured at these redshifts cannot be solely
attributed to the lower gas content in the galaxy which is much more moderate,
and we need to consider whether the star formation efficiency of the gas is
reduced compared with the HD and MHD runs. To explore the physical conditions
of the ISM gas in nut, Fig. 5 shows temperature-density phase diagrams for the
three runs used in this study. Phase diagrams based on a single simulation
snapshot can be hard to interpret as they are subject to the stochasticity
driven by processes such as feedback events and mergers. To avoid this issue,
the phase diagrams in Fig. 5 are averaged over multiple snapshots within the
galaxy’s dynamical time at $z\sim 2.5$, weighted by their difference in time
to $z\sim 2.5$. We estimate the galaxy’s dynamical time $\tau_{\rm dyn}$ as
the time required for a test particle to complete one full orbit at
$r=0.2R_{\rm vir,DM}$ with circular velocity:

$\displaystyle v_{\rm circ}(r)=\sqrt{\frac{GM(<r)}{r}}\,,$ (8)

where $M(<r)=M_{\rm b}+M_{\rm DM}$ is the combined baryon and DM mass within
that sphere. Additionally, to ensure the comparison between runs occurs over
the same time span, we use the largest $\tau_{\rm dyn}$ of the three runs
compared, which for $z\sim 2.5$ corresponds to the CRMHD simulation and has a
value of $\delta t=0.556$ Gyr. Therefore, our measurements have a central
redshift $\langle z\rangle$ and a physical time range $\delta t$.

For the HD run, the left panel of Fig. 5 shows that the low density gas
([print-unity-mantissa = false]e-26e-24 ^-3) in the galactic region is
dominated by gas with temperature $T\sim 10^{4}$ K (where hydrogen
recombination and resonance lines of atomic hydrogen dominate the cooling of
gas) and above, a typical signature of SNe shock-heating. Gas accumulates in
the $T\sim 10^{4}$ K region due to thermal instability (Field et al., 1969;
Wolfire et al., 1995). When this gas is perturbed by processes such as shock
compression, it collapses into cold dense clumps that populate the tip of the
phase diagram in the bottom right corner. In agreement with the results found
by Martin-Alvarez et al. (2018, 2020) for nut, the addition of magnetic fields
helps sustain a somewhat larger fraction of gas in hotter phases, as well as
allows for gas cells to reach lower densities whilst retaining their
temperature. This is reflected in the size of the ‘dark’ regions of the phase
diagram with densities below $10^{-24}\,\rm g\,cm^{-3}$, which are larger in
the MHD run as opposed to the HD run, due to the extra magnetic pressure
support. To further aid this analysis, we have separated the temperature-
density diagrams into three phases according to specific entropy $s$: cold and
dense (blue region), warm (orange region) and hot and diffuse (red region).
This separation is inspired by the work of Gent (2012). Similar amounts of gas
are found in the hot and warm regions of the HD and MHD runs, with the MHD run
having a slightly higher gas mass in the warm phase (from 49% in HD to 55% in
MHD). The hot phase only contributes 3% to the total gas mass in both HD and
MHD.

However, the gas phase diagram for the CRMHD run (right panel of Fig. 5) shows
marked differences compared with the non-CR models: (a) at densities
$\rho<10^{-24}\,\rm g\,cm^{-3}$, gas is no longer confined to temperatures
above $10^{4}$ K when $\rho<10^{-24}\,\rm g\,cm^{-3}$, which opens a new
region in phase space for low density cold gas; (b) the bottom right of the
warm region (i.e. $T<10^{4}$ K with $10^{-24}<\rho<10^{-22}\,\rm g\,cm^{-3}$)
has a larger scatter around the median (dashed white line) of the 2D
histogram; and (c) the cool dense region below the diagonal blue line which
separates the cold and warm phases contains less gas mass than the HD and MHD
runs. In fact, when computing the total masses in the three phases, the cold
phase mass in the CRMHD run is reduced by a factor of two compared to the HD
model as this gas now resides in the warm medium. Notably, it is in these
dense regions of the phase diagram that we expect gas to undergo a more
efficient star formation. This considerable reduction of cold gas mass
explains, at least partially, the causal link with the measured decrease of
stellar mass in this simulation as inferred from Fig. 3. Note that we have not
yet quantified the turbulent or magnetic support of the gas, and how these
affect the efficiency of star formation. Indeed, in combination with the
thermal pressure, the additional CR pressure is able to support gas against
further collapse, and the build up of a pressure gradient allows for the warm
gas reservoir above the disk to be kept there for longer (e.g. Dashyan &
Dubois, 2020; Chan et al., 2022; Farcy et al., 2022; Armillotta et al., 2022).
This depletion of cold gas, now transferred to the warm phase in the presence
of CRs is in good agreement with other works (e.g. Dubois et al., 2019; Nuñez-
Castiñeyra et al., 2022). By considering the medians of the 2D phase diagrams,
two main features are readily visible: (i) the larger contribution to the mass
budget by gas between [print-unity-mantissa = false]e-24e-22 ^-3 in the CRMHD
run creates a more pronounced concave inflection point towards lower $T$ at
$\rho\sim\qty[print-unity-mantissa=false]{e-23}{\,\centi^{-3}}$; (ii) gas mass
accumulates more tightly around the thermal instability shoulder at densities
below e-24 ^-3, driving the median line closer to this feature of the phase
diagram as the relative amount of gas above $T\geq\qty[print-unity-
mantissa=false]{e5}{}$ is greatly reduced, particularly at higher densities.

###  3.3 Cosmic rays reshape how supernovae interact with the galactic region

####  3.3.1 Distribution of supernova sites and feedback effectiveness

Figure 6: Distribution of gas densities in which SN events take place in the
redshift range $z=6.0-3.0$ (top panel) and $z=3.0-2.0$ (bottom panel). Data
corresponding to the CRMHD simulation is shown in green, the MHD simulation in
purple and the HD simulation in blue. Across all redshifts, the majority of SN
events in the HD and MHD simulations are found at high densities ([print-
unity-mantissa = false]e-22e-21 ^-3), in which SNe explosions suffer from
strong radiative cooling. Contrarily to HD and MHD, in the CRMHD simulation at
early times (top panel) the distribution is strongly bimodal, with a peak at
low densities $\sim\qty[print-unity-mantissa=false]{e-27}{\,\centi^{-3}}$ and
a peak at $\sim\qty[print-unity-mantissa=false]{e-23}{\,\centi^{-3}}$. SN
events taking place at lower densities during this period are more effective
in suppressing star formation. This bimodality for CRMHD is still present in
the lower redshift range (bottom panel), but the contribution of the low
density peak is significantly suppressed, with the overall distribution
resembling more the one obtained for the HD and MHD simulations.

The modifications we observed to the gas distribution across ISM phases also
have an effect on the distribution of SN explosions across the multi-phase
ISM. These local properties will affect the efficiency of SN feedback on self-
regulating star formation and modifying the ISM properties. We explore this
issue by recording the local ISM conditions for each SN event. In Fig. 6 we
plot the distribution of SN site gas density for two redshift ranges: between
$3<z<6$ (top panel), i.e. the period of largest $M_{*}$ difference between
CRMHD and HD/MHD, and between $2<z<3$ (bottom panel), which samples the period
of fast growth of $M_{*}$ in CRMHD. We find for the HD simulation that (as in
previous works Kimm et al., 2015) these PDFs are intrinsically bimodal with
the two peaks at low and high densities, characterising the hot and cold
phases of the ISM (see e.g. Gent, 2012). However, the high density peak
dominates the overall distribution, which at higher redshifts (top panel)
resembles more a log-normal distribution of the cold ISM. When SNe explode in
the HD case, young star particles have not had the chance to leave the dense
regions of the galaxy in which they were born. This means the $10^{51}$ erg
from a single SN888At high densities the Sedov-Taylor phase will hardly be
captured by our $\sim 20$ pc resolution, so the mechanical feedback method
from Kimm et al. (2017) will mainly inject the momentum reached at the end of
the Sedov-Taylor phase. will be deposited into a larger gas mass per cell,
hence reducing the maximum launch velocity and making it harder for the
outflowing gas to escape the gravitational potential.

Including magnetic fields in the MHD simulation moves the low density peak of
the distribution to lower densities ($\sim\qty[print-unity-
mantissa=false]{e-27}{\,\centi^{-3}}$) and produces a more bimodal
distribution at higher redshifts. As already implied by the phase diagrams of
MHD (see Fig. 5), the additional magnetic support alters the gas phase
distribution inside the galaxy, allowing for SNe to take place in lower
density gas. Furthermore, we also attribute this peak at lower densities to
the increased clustering of star-forming regions (e.g. Hennebelle & Iffrig,
2014; Hennebelle et al., 2022), as magnetic fields reduce the cold gas
fragmentation, making feedback events more correlated to each other and hence
promoting the formation of super-bubbles.

The presence of this bimodality in SN sites is strongly altered for the CRMHD
simulation compared to no-CR runs. Although the effects seen in the MHD
simulation, which arise from the support of magnetic fields, should also be in
place in the CRMHD simulation, we have already seen (in Fig. 5) that the
presence of CR pressure and heating mechanisms dominates over the effect of
magnetic fields. The ISM in the CRMHD run is strongly depleted of high density
cold gas, with the high density peak of the bimodal distribution now at
$\sim\qty[print-unity-mantissa=false]{e-23}{\,\centi^{-3}}$ in the CRMHD run
instead of $\sim\qty[print-unity-mantissa=false]{e-21}{\,\centi^{-3}}$ in the
HD and MHD runs. The regions of star formation (black contours) in Fig. 5 are
pushed to lower densities, and the warm phase dominates the gas mass budget in
the galactic region. Therefore, it is more common for SN explosions to take
place at lower densities ($\leq 10^{-22}\,\rm g\,cm^{-3}$), with a large
number of SNe exploding at $\rho\sim 10^{-27}\,\rm g\,cm^{-3}$. We attribute
this to star formation happening in a less clustered form throughout the
galactic region, with the consequent feedback events being able to sample
lower density regions shaped by CR pressure and previous SN-driven bubbles. In
addition, SN remnants evolve for longer and maintain higher temperatures past
the canonical evolutionary sequence (e.g. Naab & Ostriker, 2017) when the
effects of CR pressure are taken into account (see Rodríguez Montero et al.
2022 for a parameter space study of CR feedback effects on the evolution of
individual SN using 3D MHD simulations). Furthermore, similarly to the
behaviour of fully thermal SN remnants (Thornton et al., 1998), when CRs are
accounted for in the evolution of SN explosions in low density environments,
they can reach a larger value of deposited momentum compared to high density
environments (Rodríguez Montero et al., 2022).

The bottom panel of Fig. 6 shows the PDF of SN site densities instead between
the redshifts of 3 and 2. In this lower redshift range we find that the CRMHD
PDF begins to converge to the ones obtained for the HD and MHD runs, although
maintaining a higher fraction of SNe happening at lower density (i.e.
$\sim\qty[print-unity-mantissa=false]{e-27}{\,\centi^{-3}}$). nut transitions
into a rotationally supported (i.e. ratio of rotational velocity to dispersion
velocity larger than unity) disk at $z\sim 3.5$ (Martin-Alvarez et al., 2020)
in HD and MHD simulations, while we found that for our CRMHD run this
transition occurs closer to $z\sim 2.3$. The nut galaxy is dominated by the
gas velocity dispersion between $z=6-3$, which means that CRs released by SN
events within the galactic region can easily diffuse across the ISM and affect
additional star-forming regions. However, in the lower redshift range of the
bottom panel we are instead capturing the morphological transformation into a
dynamically cold disk. With a disk morphology in place, CRs will be more prone
to escape perpendicular to the disk, having a less drastic effect in the cold
phase of the ISM at lower redshifts. We predict a dual behaviour depending on
the morphology of the galaxy, which we will explore in future work.

####  3.3.2 Influence of cosmic rays in the launching of outflows

The CRMHD simulation has both lower $M_{\rm gas}$ and $M_{*}$ at high redshift
than the HD and MHD models. This must be, at least partially, due to
differences in the outflow and inflow rates of these simulations. In this work
we focus on the analysis of outflowing gas. We measure instantaneous outflow
rates $\dot{M}_{\rm outflow}$ as

$\displaystyle\dot{M}_{\rm outflow}(r)=\frac{1}{\delta r_{\rm
shell}}\sum_{i}\rho_{i}v_{r,i}(\Delta x_{i})^{3}\,,$ (9)

where we sum over all cells with positive radial velocities $v_{r,i}$ in a
thin spherical shell of width $\delta r_{\rm shell}$ at a distance $r$ from
the galactic centre. The width of the shell is taken to be $0.01\,R_{\rm
vir,DM}$, which provides a consistent measure across simulations and cosmic
time. Each cell is characterised by its gas density $\rho_{i}$, the radial
velocity $v_{r,i}$ and the cell width $\Delta x_{i}$. Previous studies with
the nut galaxy (Tillson et al., 2015) found the contribution to the inflowing
gas from satellites to be small ($\sim 1\%$) compared to that from filaments.
Nevertheless, for the analysis of the outflowing gas, satellites that cross
our thin shell may contain gas that could appear to be flowing radially
outward in the frame of the central galaxy (either due to satellite rotation,
or because a given satellite crosses the measuring shell on its way back to
the outer halo regions after its apocentric passage). To avoid these
subtleties, we subtract all contributions in the thin shell from passing
substructures. We thus remove gas located within the tidal radius of DM sub-
halos containing more than 1000 particles. The tidal radius is estimated as
the Jacobi radius with an isothermal sphere approximation for DM halos given
by Eq. 8.107 of Binney & Tremaine (2008). Finally, we divide outflowing gas
into separate temperature phases as follows:

  * •

Hot gas: $T>10^{5}$ K,

  * •

Warm ionised gas: $10^{5}>T>9000$ K,

  * •

Warm neutral gas: $9000>T>1000$ K and

  * •

Cold gas: $T<1000$ K .

Figure 7: Top panel: Outflow rate measured in a thin spherical shell at
$0.2\,R_{\rm vir,DM}$ for the HD, MHD and CRMHD runs. The histograms show the
instantaneous mass-flow rate through the thin shell, while the solid lines are
the running median of these measurements using a filter of size 200 Myr.
Middle panel: Star formation rate (SFR) averaged over time bins of 100 Myr.
Bottom panel: Mass-loading factor $\eta$ computed using the running median
results of $\dot{M}$ and SFR shown in the top panels. The simulation with CRs
has large outflow rates, $\times 10-100$ higher than the HD and MHD
simulations, and maintains them across all cosmic time explored. Given the
large suppression of star formation for $z\geq 3$, the CRMHD simulation has
outflow mass-loading consistently exceeding unity, only reached by the HD and
MHD simulations during the merger-induced starburst at $z\sim 3-2$.

In Fig. 7 we show the results of this outflow analysis in a thin spherical
shell at $0.2\,R_{\rm vir,DM}$ for the three simulations across $\sim 4.5\,\rm
Gyr$ of evolution, where $\dot{M}_{\rm outflow}$ (top panel), star formation
rate (SFR) averaged over $100\,\rm Myr$ (middle panel), and the mass loading
factor $\eta=\dot{M}_{\rm outflow}/{\rm SFR}$ (bottom panel) is shown. Solid
lines are the result of a running median smoothing with a size of 200 Myr.
Note that we do not restrict our analysis to particular assumptions about the
possible time delay between a star formation event in the galaxy and the
consequent outflow at $0.2\,R_{\rm vir,DM}$. We distinguish three distinct
time phases when comparing the HD, MHD and CRMHD simulations: i) between
$z\sim 13$ and $z\sim 7$ CRMHD reaches an outflow rate of
$2-4\,\rm\text{M}_{\odot}\,yr^{-1}$, whereas the HD and MHD runs are closer to
$0.1-0.4\,\rm\text{M}_{\odot}\,yr^{-1}$; ii) between $z\sim 7$ and $z\sim 3$
the outflow rates of the HD and MHD runs reach
$1\,\rm\text{M}_{\odot}\,yr^{-1}$ (with the MHD run having overall larger
values than HD), during times of strong starbursts. However, given their
average SFR of $10\,\rm\text{M}_{\odot}\,yr^{-1}$, their mass loading factors
remain well below unity, while the CRMHD simulation consistently exceeds
$\eta\sim 1$ during starbursts; and iii) between $z\sim 3$ and $z\sim 1.5$ the
SFR of the three runs peaks at $z\sim 2.5$ at a value of
$10\,\rm\text{M}_{\odot}\,yr^{-1}$ (and slowly drops to
$2-4\,\rm\text{M}_{\odot}\,yr^{-1}$). After this time, only CRMHD maintains a
steady outflow rate of $8-9\,\rm\text{M}_{\odot}\,yr^{-1}$ while MHD drops to
$2\,\rm\text{M}_{\odot}\,yr^{-1}$ and HD is kept at
$2-4\,\rm\text{M}_{\odot}\,yr^{-1}$, with a jump to
$8\,\rm\text{M}_{\odot}\,yr^{-1}$ around $z\sim 2$. It is also important to
note that in comparison with the HD and MHD runs, the CRMHD runs experiences a
fast rise in SFR of 1 order of magnitude between the redshifts of 3 and 2.3.
As mentioned before, during this period nut experiences the merger with a
large companion, causing a large increase in the cold gas mass of the galaxy
(see Fig. 3) and a burst in star formation. We find that this period is
connected with a large increase of inflow rate to nut in the CRMHD run
(Rodríguez Montero et al., in prep.), caused by the cold, high density
extraplanar gas supported by CR pressure to become unstable during the merger
interaction. Overall, the cosmic evolution of $\dot{M}_{\rm outflow}$ appears
more steady than the SFR for the CRMHD run. This suggests that while SN
feedback is rarely efficient in the HD and MHD runs with large amount of gas
consumed into stars, in the CRMHD run large quantities of gas are ejected from
the galaxy and gas-to-stars consumption rate is much lower.

Figure 8: Histogram of mass outflow rates as a function of gas radial velocity
measured in a thin shell at $0.2\,R_{\rm vir,DM}$ and $\langle z\rangle=3.94$
(top panel) and $\langle z\rangle=2.46$ (bottom panel). We compare the HD
(blue), MHD (purple), and CRMHD (green) simulations, stacking the histogram
values in snapshots spanning a $t_{\rm dyn}$ range, using the method detailed
in Section 3.2. Vertical dotted lines show the stack-weighted averages $v_{\rm
esc}$ for the redshift quoted, with the corresponding standard deviations
denoted by the shaded regions. Values for the total mass outflow rate above
$v_{\rm esc}$ ($\dot{M}_{\rm esc}$) are included in each panel for each
simulation. The peak of the outflowing gas distribution is both more prominent
and displaced to lower radial velocities in the CRMHD simulation in both
redshift ranges compared to HD and MHD. At high redshift, the shallower
gravitational potential in the CRMHD simulation allows a much larger fraction
of the high velocity tail to exceed the local $v_{\rm esc}$.

In order to further understand the underlying differences between outflows in
the non-CR runs and CRMHD, we compare their outflow velocity distributions
during two distinct phases of the SFH: 1) a peak in SFR present in all
simulations at $z\sim 4$; and 2) a period of similar SFRs across all runs
($\sim 8-10\,\rm\text{M}_{\odot}\,yr^{-1}$) occurring approximately at $z\sim
2.5$. The latter is characterised by a merger with a small companion. We
compute stacked mass-flow rate histograms following the same method detailed
in Section 3.2 for weighted stacking of PDFs, employing snapshots within a
dynamical time of the galaxy. These histograms show the total $\dot{M}$
separated into radial velocity bins for the same thin shell employed in Fig. 7
at $0.2\,R_{\rm vir,DM}$. We show the HD, MHD and CRMHD simulations in blue,
purple and green, respectively. The top panel of Fig. 8 shows that for
$z=\langle 3.94\rangle$ the outflowing gas carrying the majority of the mass
outside of the galaxy is found at lower radial velocities in the CRMHD run
than in HD and MHD. While the HD distribution is dominated by gas between
60$\,\rm km\,s^{-1}$ and 120$\,\rm km\,s^{-1}$, it rapidly drops towards lower
velocities, with a small fraction of the total outflow rate
($0.05\,\rm\text{M}_{\odot}\,yr^{-1}$) exceeding the $v_{\rm esc}$ of the
galaxy (given by the vertical dotted blue line). This escape velocity is
computed for each time $t$ by counting the total mass enclosed within
$0.2\,R_{\rm vir,DM}$. The MHD simulation instead shows a tail that extends to
higher velocities, contributing 10 times more to the mass outflow rate above
$v_{\rm esc}$ than HD. However, despite this higher efficiency driving fast
outflows, the majority of the outflowing mass between 60$\,\rm km\,s^{-1}$ and
120$\,\rm km\,s^{-1}$ remains the same, and the mass above $v_{\rm esc}$ is 3
times lower than the CRMHD. This is due to the fact that $v_{\rm esc}$ is half
the value for CRMHD than for HD/MHD, a signature of a shallower potential
well. This is due to the CR feedback significantly affecting the growth of the
central stellar mass of nut at high redshift. Therefore, the larger outflow
rates of CRMHD at $z\sim 4$ compared to HD and MHD can be explained by a
larger fraction (almost 1 dex compared to HD and MHD) of ‘low-velocity’ gas
and an easier escape of fast outflows in a less deep gravitational potential.
To further explore these trends, we also analyse the outflow velocity
distribution at $z\sim 2.5$ (see bottom panel of Fig. 8). Although now the
fast growth of nut in the CRMHD simulation between $z=3$ and $z=2$ somewhat
bridges the gap in $v_{\rm esc}$, this simulation still has 6 times as much
gas outflowing above $v_{\rm esc}$ than the HD and MHD runs. This is not
surprising given that during this period $\dot{M}_{\rm outflow}$ in CRMHD is
more than an order of magnitude higher than HD and MHD. Furthermore, these
plots also suggest that the majority of the outflowing gas at $0.2\,R_{\rm
vir,DM}$ lacks the kinetic energy necessary to escape the local gravitational
potential in all simulations. However, it remains an open question whether
this gas needs to have the energy necessary to escape the halo in order to
reduce the reservoir of gas available for star formation or if it is enough to
prevent its re-accretion over sufficiently long timescales (see e.g.
Somerville & Davé, 2015, for preventive feedback in the context of AGNs). As
far as our current analysis goes, the CRMHD model is able to reach higher mass
loadings than the no-CR runs. These outflows are mainly comprised of a slow
moving wind that, without any further acceleration beyond $0.2R_{\rm vir,DM}$,
will not be able to escape the halo and will eventually fall back into the
galaxy if undisturbed by future feedback events.

####  3.3.3 Dynamics of multi-phase outflows

Figure 9: Outflow budget and kinematics of hot (top left pair), warm ionised
(top right pair), warm neutral (bottom left pair) and cold (bottom right pair)
phases for the HD (blue), MHD (purple) and CRMHD (green) runs measured at
$0.2\,R_{\rm vir,DM}$. Each pair of panels show the instantaneous outflow rate
(top) and the $\dot{M}_{\rm outflow}$-weighted radial velocity of outflowing
gas (bottom). The $\dot{M}_{\rm outflow}$ panels show the same quantities as
Fig. 7: low opacity histogram is the instantaneous value, while the solid
lines show the running median with a filter of $200$ Myr. In the $v_{r}$
panels, the escape velocity of each model at the shell radius is showed with
dashed lines, for ease of comparison to the median $v_{\rm outflow}$. Shaded
regions represent the second and fourth quartiles of the underlying weighted
distribution. Hot gas attains the largest radial velocities of all phases,
with an almost constant value of $\sim 100\,\rm km\,s^{-1}$ for $z\geq 3$ and
$\sim 70\,\rm km\,s^{-1}$ for $z\lesssim 3$, dominating the outflow rate in HD
and MHD simulations. The hot phase in CRMHD reaches the largest outflow
velocities ($\sim 150-200\,\rm km\,s^{-1}$ for $z\geq 2.5$ and $\sim 100\,\rm
km\,s^{-1}$ for $z\lesssim 2.5$), with a non-negligible fraction of the
outflowing gas mass exceeding $v_{\rm esc}$ consistently up to $z\geq 2.5$.
The dominant mass-flow rate is coming from the warm ionised phase in the
CRMHD, while above $z\sim 5$, warm neutral and cold phases largely contribute
to the total outflow rate. Warm neutral and cold phases are scarce in the HD
and MHD simulations, usually identified as low velocity peaks which can
matched with crossing gas features of orbiting satellites. However, in the
CRMHD run these outflows while slower than the hot gas can still reach $\sim
100\,\rm km\,s^{-1}$ at lower redshifts.

The higher outflow rates measured for the CRMHD simulation despite a smaller
SFR are a consequence of the energy deposited by SNe being more efficient at
expelling gas from the galactic disk. To provide evidence for this higher
efficiency causing these higher mass-loading results in CRMHD, Fig. 9 shows
the $\dot{M}_{\rm outflow}$ and $\dot{M}_{\rm outflow}$-weighted radial
velocities of outflowing gas $v_{\rm outflow}$ in the HD (blue), MHD (purple)
and CRMHD (green) simulations, separated in the different phases we explore.
HD and MHD outflows are dominated by the hot and ionised phases during the
majority of the evolution of nut. Outflowing gas in the warm neutral and cold
phases in these runs are identified with sudden, low velocity peaks in the
velocity plots, which we found to be related to clumps of gas in eccentric
orbits around nut, rather than genuine outflows. Instead, the CRMHD run shows
that the hot phase never dominates the outflowing mass rate. The warm ionised
phase dominates ($1-10\,\rm\text{M}_{\odot}\,yr^{-1}$) the outflowing gas mass
up until $z\sim 4$, with considerable contributions of warm neutral and cold
gas at higher redshifts. We find that multi-phase outflows ($\sim$1:2:4:2
ratios for cold, warm neutral, warm ionised and hot phases, respectively) are
only found consistently in the CRMHD simulation, and at high redshift ($z\geq
4$). The HD and MHD runs fail to produce such gas mass distributions across
phases at all redshifts presented here.

We compare the radial velocities with the escape velocity at $0.2\,R_{\rm
vir,DM}$. To understand whether the outflowing gas is gravitationally bound,
we compare these velocities with the local escape velocity shown in the four
panels with dashed lines. While a radial velocity higher than $v_{\rm esc}$ is
not a sufficient condition for gas to completely escape the gravitational
potential of a halo, it is a useful proxy commonly used in observations and
theoretical studies. The radial velocities of the hot phase (top left panel)
range from $\sim 100\,\rm km\,s^{-1}$ for $z\geq 3$ to $\sim 80\,\rm
km\,s^{-1}$ at lower redshifts in the no-CR simulations. For the case of
CRMHD, the velocities are consistently higher, being in the range $\sim
100-200\,\rm km\,s^{-1}$. These outflow-weighted velocities are in agreement
with observations of the warm ionised gas in local dwarf galaxies by McQuinn
et al. (2019) and Marasco et al. (2023), while falling short of the $\sim
400-600\,\rm km\,s^{-1}$ value reported in Concas et al. (2022) for their
sample of galaxies between $1.2<z<2.6$. However, the CRMHD simulation is the
only one consistently able to reach outflow velocities above the escape
velocity within the galactic region, up until $z\sim 2.5$, when its
distribution becomes similar to the non-CR cases. This behaviour is in
agreement with the histograms shown in Fig. 8.

Outflows in the warm neutral and cold phases are rarely found in the HD and
MHD runs (except during the merger taking place at $z\sim 3-2$), and only
appear as rare and sharp peaks in the $v_{r}$ plots. By examining projections
of the galactic region of nut at the redshifts corresponding to these peaks,
we have found that for HD and MHD all cold and warm neutral outflows are
linked to extended tidal features of satellites, and are not spatially related
to the rest of hot and warm ionised outflowing gas. In fact, CRMHD is the only
run that shows warm neutral and cold outflows across the majority of
snapshots. These outflows are also spatially correlated with the hotter phases
of the outflowing gas and are not dynamically linked with tidally stripped
satellites (see Section 3.3.4 and Fig. 11 for further details). Contrarily to
the hot and warm ionised phases, the warm neutral and cold gas appear to show
an evolution in outflow velocity, increasing from $\sim 50\,\rm km\,s^{-1}$ at
high redshift to a maximum of $\sim 100\,\rm km\,s^{-1}$ at $z\sim 2.5-2$.
Similarly to the warm ionised phase, the median $v_{r}$ of these two phases is
below $v_{\rm esc}$ (except for peaks at $\sim 2.2$ and $\sim 2.5\,\rm Gyr$ in
the cold phase). Overall, we find that these colder outflows struggle to reach
the escape velocities at $0.2\,R_{\rm vir,DM}$, even in the CRMHD run and
despite its shallower gravitational potential.

####  3.3.4 Cosmic ray-driven re-acceleration of outflows in the CGM

Figure 10: Pressure gradient support against gravity for outflowing gas
separated into different temperature phases for HD (blue), MHD (purple) and
CRMHD (green) simulations. The dividing line between the shaded and unshaded
region indicates marginal support. In the shaded region, gravity dominates the
dynamics, while in the unshaded region, gas is accelerated outwards by the
pressure gradients. Outflowing gas is measured in a thin shell at $0.2\,R_{\rm
vir,DM}$ and divided into the phases described in Section 3.3.2 according to
temperature. Shaded, coloured bands represent the second and fourth quartiles
of the underlying distribution. The left column corresponds to the hot gas,
with each row corresponding to the HD (top), MHD (middle), and CRMHD (bottom)
simulations, respectively. While below $z\sim 3$ the thermal pressure
struggles to accelerate gas in the HD and MHD simulations, CR pressure
consistently supports and accelerates the hot gas across all cosmic time. Only
CR pressure can mildly support the warm ionised phase between $3<z<6$, but it
supports and accelerates warm neutral and cold gas for $z\leq 5$.

The CRMHD simulation displays considerably higher outflow rates at
$0.2\,R_{\rm vir,DM}$ than the HD and MHD runs. However, the previous analysis
has shown that gas usually lacks the kinetic energy required to fully escape
the gravitational potential of the galaxy. In addition, outflowing gas is
further affected by the interaction with the CGM. Therefore, in this section,
we investigate how such an interaction further impacts the gas dynamics at
$0.2\,R_{\rm vir,DM}$ in Fig. 10. This is carried out by computing the mass-
flow-rate-weighted median of the radial pressure support for the outflowing
gas in the thin shell at $0.2\,R_{\rm vir,DM}$. As per the hydrostatic
equilibrium equation, the negative pressure gradient of the gas needs to
balance the inward pull of gravity for a parcel of gas not to inflow. The
gravitationally-dominated regime of which this equilibrium situation is a
limiting case, is indicated by the grey band below a value of 1. For a
pressured fluid $x$, this corresponds to
$-\hat{r}\cdot\vec{\nabla}P_{x}=\hat{r}\cdot(\rho\vec{g})$, where $\hat{r}$ is
the radial unit vector, $\vec{g}$ is the local gravitational acceleration and
$x$ is the fluid pressure component (either thermal or CR) considered. Cells
with radial support below unity experience an overall deceleration whereas
those with a value $>1$ experience an outward acceleration. We represent the
contributions from thermal pressure, $P_{\rm ther}$, as solid circle symbols
coloured in red linked by dashed lines, and CR pressure, $P_{\rm CR}$, as
green circles connected by solid lines.

We begin by comparing the support for the hot phase, where we separate for
clarity the HD (top left), MHD (middle left) and CRMHD (bottom left)
simulations. We found that the hot and warm ionised phases show a large
variance for these quantities, so we include the second and fourth quartiles
of the underlying distribution as a shaded band in order to exemplify this.
The three runs display a large amount of thermal pressure support against
gravity for $z\geq 3$. At lower redshifts this pressure support ratio rarely
exceeds unity (although it is more common for the MHD case), and has a more
bursty behaviour. This change in the support of hot outflowing gas in all
simulations is attributed to the growth of the virial shock at lower
redshifts, as the cooling rate becomes smaller than the post-shock compression
rate (Dekel & Birnboim, 2006) and cold accretion becomes less important
(Powell et al., 2011). A hotter CGM means a larger thermal pressure confining
gas to the galactic disk, and buoyancy forces present at high redshift becomes
less important. Although, this behaviour is observed for all runs, the CRMHD
run features a larger frequency of radial support peaks exceeding the force
balance threshold, despite the higher density of the CGM in the CRMHD
simulation compared with the non-CR simulations (see Section 3.1). We
attribute this to a colder CGM with a lower fractional contribution from
thermal support due to the additional support of CR pressure (see Section
3.1). We do not include the contribution of magnetic support, as we find it to
be negligible for outflows at all times for the MHD and CRMHD
simulations999Nevertheless, this should be considered also in the context of
magnetic field injection by SN remnants (e.g. Martin-Alvarez et al., 2021),
which we do not explore in this study.. Thermal pressure dominates the hot
phase in CRMHD down to $z\sim 3$, with both CRs and thermal pressures
separately capable of supporting a considerable fraction of the outflowing gas
against gravity. The CR pressure support of the hot component evolves only
mildly with redshift, eventually dominating for $z<3$. This change in
behaviour coincides with the drop in hot gas radial velocity after $z\sim 3$
shown for CRMHD in Fig. 9. While hot fast outflows in the HD and MHD
simulations find it harder to travel through a the CGM at lower redshifts, in
the CRMHD the additional support of CRs injected by SNe to the wind bubble
further accelerates them when their thermal buoyancy stalls.

Figure 11: Edge-on view of nut in the CRMHD run during a strong feedback event
at $z=2.41$. Top panels: density-weighted projections of a $20\,\rm kpc$ on a
side cube centred on nut showing the gas radial velocity (left panel) and the
gas temperature (right panel). In the left panel, inflowing gas is shown in
blue and outflowing gas in red. The temperature map indicates that
considerable hot outflows are ejected perpendicular to the cold disk, while
warm gas at $\sim 10^{4.5}\,\rm K$ surrounds the galaxy. Black (green) and red
(blue) contours are overlaid on the radial velocity (and temperature)
projection to show gas density contours at $10^{-24}\,\rm g\,cm^{-3}$ and
$10^{-22}\,\rm g\,cm^{-3}$, respectively. The purple dashed circle shows the
$0.2\,R_{\rm vir,DM}$ shell used to measure gas flows close to the galactic
region. Bottom panels: the same projected region now separated into four
outflow gas phases (see text for detail). Each phase is plotted in two panels,
corresponding to thermal (i.e. $-\nabla_{r}P_{\rm ther}/(\rho g_{r})$) and CR
(i.e. $-\nabla_{r}P_{\rm CR}/(\rho g_{r})$) support of the gas against
gravity, with red (outwards acceleration) or blue (inwards acceleration)
colours indicating the dominant pressure gradient. From top left to bottom
right, the temperature phases displayed are: hot (two top leftmost panels),
warm ionised (two top rightmost panels), warm neutral (two bottom leftmost
panels), and cold (two bottom rightmost panels). To help visual analysis, we
include the density contours (red and blue) of the temperature panel in these
panels. The warm neutral and cold panels also include contour plots showing
the regions dominated by the hot (blue contours) outflows. Hot outflows are
clearly biconical with large gradients of both thermal and CR pressure. Both
feature alternative positive and negative gradients, a tell-tale signature of
consecutive outflow waves. The warm ionised gas is the most volume filling
phase, with a much more isotropic structure around nut, and dominated by a
uniform $-\nabla_{r}P_{\rm CR}/(\rho g_{r})$ profile. Warm neutral and cold
outflows present a more clumpy and filamentary structure which quickly fades
off at larger radii.

We compare in the right panels of Fig. 10 the pressure support against gravity
in the warm ionised, warm neutral and cold phases, now for the three
simulations combined. As warm and cold outflows are scarce in the HD and MHD
runs, we limit our analysis for these phases to the CRMHD run. These confirm
that the thermal pressure support (as well as the magnetic, which we do not
include) in the HD and MHD runs is negligible. This further supports the idea
that for the runs without CRs, outflowing gas colder than $10^{5}$ K does not
experience further acceleration at $0.2\,R_{\rm vir,DM}$. Due to outflow
velocities lower than $v_{\rm esc}$, outflows likely rain back down onto the
galaxy as a galactic fountain. The CRMHD run features warm ionised outflows
that are marginally supported against gravity, following an overall trend
similar to that of the hot gas, but now with CR pressure dominating over
thermal pressure at all times. In fact, CR pressure support dominates over
thermal pressure for the three phases in the right column of Fig. 10. The warm
neutral and the cold components of the outflow show large accelerations due to
CR pressure, exceeding by up to 7 times the gravitational acceleration. These
phases display a large variability, also found in the analysis of $v_{\rm r}$
(see Fig. 9). This can be attributed to a combination of outflow ‘burstiness’
and a lower volume-filling of these phases. Therefore, whenever outflows
entraining warm neutral and cold gas are present, they are usually subject to
large CR pressure gradients that not only balance gravity but provide
additional acceleration. This means that, although the outflow velocities of
these two phases at these radii are not enough to escape the halo, CR pressure
gradients can drastically change their dynamics as winds escape from the
galaxy.

This section has described how CRs affect outflows, providing further insight
on the consistently higher outflow rates that the CRMHD run exhibits when
compared to the HD and MHD runs. Here we perform a more detailed analysis of a
particular feedback event in the CRMHD simulation to understand how CRs are
locally driving the observed differences. In Fig. 11 we show edge-on density-
weighted projections of nut during a feedback event at $z=2.41$ ($t\sim
2.9\,\rm Gyr$). This time was selected because of the presence of a
corresponding peak in the outflow rate (i.e. $\sim
30\,\rm\text{M}_{\odot}\,yr^{-1}$, see Fig. 7) and the CR-supported warm
neutral and cold outflows. The top panels show the radial velocity (left) and
gas temperature (right) maps for a region $20$ kpc on a side, encompassing the
$0.2\,R_{\rm vir,DM}$ thin shell (shown as a dashed purple circle). We overlay
on top of the left (right) map contours for the gas density showing [print-
unity-mantissa = false]e-24 ^-3 in black (red) and at [print-unity-mantissa =
false]e-22 ^-3 in red (blue). These panels uphold the canonical understanding
of wind launching in disk galaxies: at the centre lies a cold and dense disk
from which high velocity, hot outflows carve biconical perpendicular channels
through the warm ionised CGM. This particular outflow is preceded by previous
plumes of hot gas visible at larger radii. There are also regions of
outflowing gas that do not have a significantly higher temperature than the
surrounding CGM (see outflowing region in the top left corner of the $v_{\rm
r}$ map).

In keeping with our previous analysis, we separate the outflowing gas into
distinct temperature phases in the 8 lower panels. For each of the phases, we
plot radial support against gravity, i.e. $-\nabla_{r}P_{x}=\rho g_{r}$, for
thermal ($P_{\rm ther}$) and CR ($P_{\rm CR}$) pressures. To help the
analysis, we keep the density contours (green and blue lines) of the
temperature projection in these panels. We first examine the hot outflows (top
left panels), which display in both thermal and CR pressure support maps a
clear biconical morphology with opening angle of $\sim 30^{\circ}$. Dark red
and blue colours represent high pressure acceleration anti-parallel and
parallel to the gravitational acceleration, respectively. The intertwined
structure of positive and negative gradients in both thermal and CR pressure
is attributed to the succession of multiple outflow bubbles, ejected from the
galactic disk at different times. As such, the front of a previous bubble
appears as a positive gradient in the radial direction. These features are
also present in the pressure support maps of the HD and MHD simulations (not
shown). This reaffirms the selection of the positive support part of these
outflows (see Fig. 10) as a method to identify galactic winds. In agreement
with the peak observed at $t\sim 2.9\,\rm Gyr$ in Fig. 10, the thermal and CR
supports are much larger than unity, both reaching similar magnitudes.
Nevertheless, large positive or negative thermal and CR pressure gradients do
not always spatially coincide. We surmise that such offsets could be due to CR
streaming, as on $\sim\,\rm kpc$ scales it can be an important driver of CR
transport and displace CRs with respect to the gas, while maintaining the CR
pressure gradient (e.g. Dubois et al., 2019).

The warm ionised phase of the gas is the dominant volume-filling phase of the
outflow, and has a more spherically isotropic morphology. This phase does not
display strong alternating features from consecutive wind bubbles, but rather
a smooth pressure gradient. This pressure support extends to large radii and
is dominated by CR pressure. In fact, the warm ionised, warm neutral and cold
phases exhibit the same support morphology, which suggests that they are all
sourced by a common underlying CR pressurised halo that surrounds the galaxy.
However, based on the mild support of CR pressure on the warm ionised phase in
Fig. 10, we claim that this phase will probably rain down as a galactic
fountain given the low average radial velocity (see top right panel of Fig.
9). The four bottom panels examine the warm neutral (left) and cold (right)
outflows. Gas in these phases is distributed in small structures that resemble
cooling filaments. Blue contours are included to indicate the limits of the
hot outflow cones. These colder phases appear to also be moving
perpendicularly to the disk, experiencing a uniform CR pressure gradient that
dominates their gravitational acceleration. This indicates that these cold
regions of gas are still efficiently accelerated by the surrounding hotter
phase, and hence entrained in the wind. The radial distribution of these
outflowing phases reveals that they decrease in fraction of volume occupied as
a function of radial distance, thus allowing to better showcase their clumpy
nature. These small outflowing gas structures do traverse the $0.2\,R_{\rm
vir,DM}$ shell we employ to measure outflows (illustrated by a dashed purple
circle) in a timescale shorter than our snapshot time resolution, naturally
resulting in the fluctuating behaviour of warm neutral and cold phases present
in Figs. 9 and 10.

##  4 Conclusions

In this work, we have performed a detailed study of the effects of CR feedback
on a Milky Way-like galaxy using a suite of new zoom-in cosmological
simulations. The main galaxy in our zoom region, nut, was simulated down to
$z=1.5$ with the ramses code. Using adaptive mesh refinement, we simulated the
multi-phase ISM with a spatial resolution of $\sim 23\,\rm pc$ per full cell
size. Our hydrodynamical simulation (HD) includes a model for a resolved ISM,
with a magneto-thermo-turbulent (MTT) star formation prescription and a
’mechanical’ SN feedback of a canonical strength. In our magnetohydrodynamic
simulation (MHD) we additionally followed the magnetic fields with the
constrained transport method which are amplified from a primordial seed of
$3\times 10^{-12}\,\rm G$. Our cosmic ray MHD simulation (CRMHD) further
models CRs injected by supernova (SN) events with an efficiency of 10%. In
this first work in a series of papers analysing the role of CRs in the nut
suite, we focused on the impact of our full CR physics which includes
advection, losses, anisotropic diffusion and streaming on the ISM and the
launching of outflows. Our main results can be summarised as follows.

  1. 1.

Including CRs significantly reduces the resulting stellar mass of nut compared
to the no-CR models, yielding good agreement with abundance matching models
(e.g. Behroozi et al., 2013; Moster et al., 2013) without the need for a
calibrated boost to the SN energy deposition or calibrated galactic winds.
This reduction in stellar mass is significant and it reaches approximately an
order of magnitude at high redshift $z\gtrsim 3$ before tapering to a factor
$\sim 4$ at $z\sim 1.5$ (Section 3.1 and Fig. 4).

  2. 2.

CRs pressure gradients support cold, low density gas against collapse, making
the warm phase significantly more dominant in the ISM (Section 3.2). CRs also
deplete an important fraction of the mass locked in the cold, dense phase
(i.e. $\rho\gtrsim 10^{-21}\,\rm g\,cm^{-3}$ and $T<10^{3}\,\rm K$) present in
non-CRs simulations, and hence significantly decrease the amount of gas
efficiently forming stars. In fact, galactic gas is rarely found below
$10^{-27}\,\rm g\,cm^{-3}$ and the ISM has a smoother gas distribution within
the disk.

  3. 3.

As a consequence of the thermodynamical changes to the structure of the ISM
and a CR-aided expansion of SN bubbles (e.g. Rodríguez Montero et al., 2022),
the environments where SN events take place are modified. In the CRMHD
simulation, SN explosion sites have a bi-modal distribution that both peaks at
and reaches considerably lower densities than in the HD and MHD runs (see Fig.
6). This reduces the amount of SN energy lost due to cooling, increases the
efficiency of this feedback mechanism and drives stronger outflows.

  4. 4.

The mass outflow rates in the CRMHD simulation are significantly larger at all
cosmic times with respect to the no-CR runs, specifically with a higher mass
loading factor $\eta$ especially at high redshift. In contrast to the non-CR
simulations where the dominant phase of outflowing gas mass is hot (i.e.
$T>10^{5}\,\rm K$), CRMHD outflows are dominated in mass by the warm ionised
phase (i.e. $9000<T<10^{5}\,\rm K$) and have a significant contribution from
the cold phase (see Section 3.3.2 and 3.3.3). Furthermore, hot outflows are
biconical in structure whereas the colder phases embedded in the hot outflow
are filamentary and clumpy. This is a direct prediction for upcoming JWST
observations of multi-phase galactic outflows expected soon.

  5. 5.

Outflows in the CRMHD simulation experience further acceleration due to the CR
pressure gradient, with a secondary contribution from the thermal pressure.
Outflows in the warm and cold phases have a significant support from the CRs,
while the thermal pressure support remains negligible. Therefore, outflows
experience further acceleration within the CGM of our simulated galaxy (see
Section 3.3.4 and Fig. 11).

We find CRs to have a large impact across various properties of the simulated
galaxy and its ISM. While we include multiple important CR physical components
in our simulations, the specific impact of each is not fully explored. A
review of these roles is left to the next paper in this series, which will
analyse the effects of this CR feedback on the thermodynamic state of the CGM
outside the nut galaxy (Rodríguez Montero et al. in prep.).

While simulations, such as the ones presented here, are revealing CRs to be a
central agent in galaxy formation, our CR models are far from being complete
and our results suffer from multiple caveats, both physical and numerical. For
instance, we still lack a self-consistent coupling of CR transport with the
local thermodynamic state of the gas (e.g. Armillotta et al., 2022).
Furthermore, CR physics may well be affected by more realistic and complete
gas cooling prescriptions (Nuñez-Castiñeyra et al., 2022; Katz et al., 2022)
and different refinement criteria in the CGM (Bennett & Sijacki, 2020; Rey et
al., 2023). Milky Way-like halos at high redshift can have their star
formation strongly altered by other forms of feedback such as that from
stellar winds (Agertz et al., 2021), AGN (Koudmani et al., 2022), or radiation
from young stars(Rosdahl & Teyssier, 2015; Hopkins et al., 2020; Emerick et
al., 2018). Their self-consistent interaction with CR feedback remains a topic
for further investigation (e.g. Farcy et al., 2022; Martin-Alvarez et al.,
2022).

Despite these caveats, our work highlights the fundamental role CRs likely
play in the formation and evolution of Milky Way-like galaxies. In particular,
they seem to be a central contributor to the feedback loop regulating star
formation and driving galactic outflows, providing a promising avenue to form
galaxies with more realistic physical properties without having to resort to
ad-hoc feedback model calibrations.

## Acknowledgements

FRM is supported by the Wolfson Harrison UK Research Council Physics
Scholarship. The authors thank Romain Teyssier for making the ramses code
publicly available. This project has received funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation programme (grant agreement No 693024). Simulations were performed
using the DiRAC Data Intensive service at Leicester, operated by the
University of Leicester IT Services, which forms part of the STFC DiRAC HPC
Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding
via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC
Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure.
Analysis of the simulations took place on the Oxford Astrophysics cluster
‘Glamdring’, that we warmly thank J. Patterson for running smoothly.

## Data Availability

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

## References

  * Agertz et al. (2013) Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, ApJ, 770, 25
  * Agertz et al. (2021) Agertz O., et al., 2021, MNRAS, 503, 5826
  * Armillotta et al. (2022) Armillotta L., Ostriker E. C., Jiang Y.-F., 2022, ApJ, 929, 170
  * Axford (1981) Axford W. I., 1981, Ann. NY Acad. Sci., 375, 297
  * Baldry et al. (2020) Baldry I. K., Sullivan T., Rani R., Turner S., 2020, MNRAS, 500, 1557
  * Barro et al. (2013) Barro G., et al., 2013, ApJ, 765, 104
  * Basu & Roy (2013) Basu A., Roy S., 2013, MNRAS, 433, 1675
  * Beck (2016) Beck R., 2016, A&ARv, 24, 4
  * Beckmann et al. (2022) Beckmann R. S., Dubois Y., Pellissier A., Olivares V., Polles F. L., Hahn O., Guillard P., Lehnert M. D., 2022, A&A, 665, A129
  * Behroozi et al. (2013) Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57
  * Bell (1978) Bell E. F., 1978, MNRAS, 182, 147
  * Bennett & Sijacki (2020) Bennett J. S., Sijacki D., 2020, Monthly Notices of the Royal Astronomical Society, 499, 597
  * Biermann (1950) Biermann L., 1950, Zeitschrift Naturforschung Teil A, 5, 65
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics, 2nd edn. Princeton University Press, http://press.princeton.edu/titles/8697.html
  * Blandford & Eichler (1987) Blandford R., Eichler D., 1987, PhR, 154, 1
  * Blondin et al. (1998) Blondin J. M., Wright E. B., Borkowski K. J., Reynolds S. P., 1998, ApJ, 500, 342
  * Booth et al. (2013) Booth C. M., Agertz O., Kravtsov A. V., Gnedin N. Y., 2013, ApJ, 777, L16
  * Boulares & Cox (1990) Boulares A., Cox D. P., 1990, ApJ, 365, 544
  * Bournaud et al. (2010) Bournaud F., Elmegreen B. G., Teyssier R., Block D. L., Puerari I., 2010, MNRAS, 409, 1088
  * Bournaud et al. (2013) Bournaud F., et al., 2013, ApJ, 780, 57
  * Breitschwerdt et al. (1991) Breitschwerdt D., McKenzie J. F., Voelk H. J., 1991, A&A, 245, 79
  * Brooks et al. (2007) Brooks A. M., Governato F., Booth C. M., Willman B., Gardner J. P., Wadsley J., Stinson G., Quinn T., 2007, ApJ, 655, L17
  * Brüggen & Scannapieco (2020) Brüggen M., Scannapieco E., 2020, ApJ, 905, 19
  * Buck et al. (2020) Buck T., Pfrommer C., Pakmor R., Grand R. J. J., Springel V., 2020, MNRAS, 497, 1712
  * Bustard et al. (2021) Bustard C., Zweibel E. G., Bustard C., Zweibel E. G., 2021, ApJ, 913, 106
  * Camps & Baes (2015) Camps P., Baes M., 2015, Astron. Comput., 9, 20
  * Cañas et al. (2019) Cañas R., Elahi P. J., Welker C., Lagos C. d. P., Power C., Dubois Y., Pichon C., 2019, MNRAS, 482, 2039
  * Caprioli & Spitkovsky (2014) Caprioli D., Spitkovsky A., 2014, ApJ, 794, 47
  * Chan et al. (2019) Chan T. K., Kereš D., Hopkins P. F., Quataert E., Su K.-Y., Hayward C. C., Faucher-Giguère C.-A., 2019, MNRAS, 488, 3716
  * Chan et al. (2022) Chan T. K., Kereš D., Gurvich A. B., Hopkins P. F., Trapp C., Ji S., Faucher-Giguère C.-A., 2022, MNRAS, 517, 597
  * Chisholm et al. (2016) Chisholm J., Tremonti Christy A., Leitherer C., Chen Y., 2016, MNRAS, 463, 541
  * Chisholm et al. (2018) Chisholm J., Tremonti C., Leitherer C., 2018, MNRAS, 481, 1690
  * Collins & Read (2022) Collins M. L., Read J. I., 2022, Nat. Astron., 6, 647
  * Commerçon et al. (2019) Commerçon B., Marcowith A., Dubois Y., 2019, A&A, 622, A143
  * Concas et al. (2022) Concas A., et al., 2022, MNRAS, 513, 2535
  * Crain et al. (2015) Crain R. A., et al., 2015, MNRAS, 450, 1937
  * Dashyan & Dubois (2020) Dashyan G., Dubois Y., 2020, A&A, 638, 123
  * Davé et al. (2019) Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827
  * De Lucia & Blaizot (2007) De Lucia G., Blaizot J., 2007, MNRAS, 375, 2
  * Dekel & Birnboim (2006) Dekel A., Birnboim Y., 2006, MNRAS, 368, 2
  * Dobbs et al. (2011) Dobbs C. L., Burkert A., Pringle J. E., 2011, MNRAS, 417, 1318
  * Dubois & Commerçon (2016) Dubois Y., Commerçon B., 2016, A&A, 585, 1
  * Dubois et al. (2013) Dubois Y., Gavazzi R., Peirani S., Silk J., 2013, MNRAS, 433, 3297
  * Dubois et al. (2016) Dubois Y., Peirani S., Pichon C., Devriendt J., Gavazzi R., Welker C., Volonteri M., 2016, MNRAS, 463, 3948
  * Dubois et al. (2019) Dubois Y., Commerçon B., Marcowith A., Brahimi L., 2019, A&A, 631, A121
  * Dunkley et al. (2009) Dunkley J., et al., 2009, ApJS, 180, 306
  * Durrer & Neronov (2013) Durrer R., Neronov A., 2013, A&ARv, 21, 62
  * Elahi et al. (2019) Elahi P. J., Cañas R., Poulton R. J. J., Tobar R. J., Willis J. S., Lagos C. d. P., Power C., Robotham A. S. G., 2019, PASA, 36, 26
  * Emerick et al. (2018) Emerick A., Bryan G. L., Low M.-M. M., 2018, ApJ, 865, L22
  * Ensslin (2004) Ensslin T., 2004, J. Korean Astron. Soc., 37, 439
  * Evans II (1999) Evans II N. J., 1999, ARA&A, 37, 311
  * Evans et al. (2009) Evans N. J., et al., 2009, ApJS, 181, 321
  * Everett et al. (2008) Everett J. E., Zweibel E. G., Benjamin R. A., McCammon D., Rocks L., Gallagher III J. S., 2008, ApJ, 674, 258
  * Farber et al. (2018) Farber R., Ruszkowski M., Yang H.-Y. K., Zweibel E. G., 2018, ApJ, 856, 112
  * Farcy et al. (2022) Farcy M., Rosdahl J., Dubois Y., Blaizot J., Martin-Alvarez S., 2022, MNRAS, 513, 5000
  * Farmer & Goldreich (2003) Farmer A. J., Goldreich P., 2003, ApJ, 604, 671
  * Federrath & Klessen (2012) Federrath C., Klessen R. S., 2012, ApJ, 761, 156
  * Ferland et al. (1998) Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B., Verner E. M., 1998, Publ. Astron. Soc. Pac., 110, 761
  * Field et al. (1969) Field G. B., Goldsmith D. W., Habing H. J., Field G. B., Goldsmith D. W., Habing H. J., 1969, ApJL, 155, L149
  * Fierlinger et al. (2016) Fierlinger K. M., et al., 2016, MNRAS, 456, 710
  * Fromang et al. (2006) Fromang S., Hennebelle P., Teyssier R., 2006, A&A, 457, 371
  * Geen et al. (2013) Geen S., Slyz A., Devriendt J., 2013, MNRAS, 429, 633
  * Geen et al. (2015) Geen S., Rosdahl J., Blaizot J., Devriendt J., Slyz A., 2015, MNRAS, 448, 3248
  * Gent (2012) Gent F. A., 2012, Space Sci. Rev., 166, 281
  * Girichidis et al. (2016) Girichidis P., et al., 2016, ApJ, 816, L19
  * Governato et al. (2010) Governato F., et al., 2010, Nature, 463, 203
  * Governato et al. (2012) Governato F., et al., 2012, MNRAS, 422, 1231
  * Guo & Oh (2008) Guo F., Oh S. P., 2008, MNRAS, 384, 251
  * Haardt & Madau (1996) Haardt F., Madau P., 1996, ApJ, 461, 20
  * Hanasz et al. (2013) Hanasz M., Lesch H., Naab T., Gawryszczak A., Kowalik K., Wóltański D., 2013, ApJ, 777, L38
  * Hanasz et al. (2021) Hanasz M., Strong A. W., Girichidis P., 2021, Living Rev. Sol. Phys., 7, 2
  * Heckman et al. (1995) Heckman T. M., et al., 1995, ApJ, 448, 98
  * Hennebelle & Chabrier (2011) Hennebelle P., Chabrier G., 2011, ApJL, 743
  * Hennebelle & Iffrig (2014) Hennebelle P., Iffrig O., 2014, A&A, 570, A81
  * Hennebelle et al. (2022) Hennebelle P., et al., 2022, A&A, 668, A147
  * Hopkins et al. (2011) Hopkins P. F., Quataert E., Murray N., 2011, MNRAS, 417, 950
  * 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. (2020) Hopkins P. F., et al., 2020, MNRAS, 492, 3465
  * Hopkins et al. (2021a) Hopkins P. F., Chan T. K., Ji S., Hummels C. B., Kereˇs D., Quataert E., Faucher-Giguère C. A., 2021a, MNRAS, 501, 3640
  * Hopkins et al. (2021b) Hopkins P. F., Squire J., Chan T. K., Quataert E., Ji S., Kereš D., Faucher-Giguère C.-A., 2021b, MNRAS, 501, 4184
  * Inoue et al. (2016) Inoue S., Dekel A., Mandelker N., Ceverino D., Bournaud F., Primack J., 2016, MNRAS, 456, 2052
  * Ipavich (1975) Ipavich F. M., 1975, ApJ, 196, 107
  * Ji et al. (2020) Ji S., et al., 2020, MNRAS, 496, 4221
  * Ji et al. (2021) Ji S., et al., 2021, MNRAS, 505, 259
  * Jóhannesson et al. (2019) Jóhannesson G., Porter T. A., Moskalenko I. V., 2019, ApJ, 879, 91
  * Jubelgas et al. (2008) Jubelgas M., Springel V., Enßlin T., Pfrommer C., 2008, A&A, 481, 33
  * Katz et al. (2022) Katz H., et al., 2022, arXiv e-prints, p. arXiv:2211.04626
  * Kay et al. (2002) Kay S. T., Pearce F. R., Frenk C. S., Jenkins A., 2002, MNRAS, 330, 113
  * Kennicutt (1998) Kennicutt R. C., 1998, ARA&A, 36, 189
  * Kereš et al. (2009) Kereš D., Katz N., Davé R., Fardal M., Weinberg D. H., 2009, MNRAS, 396, 2332
  * 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
  * Kimm & Cen (2014) Kimm T., Cen R., 2014, ApJ, 788, 121
  * Kimm et al. (2015) Kimm T., Cen R., Devriendt J., Dubois Y., Slyz A., 2015, MNRAS, 451, 2900
  * Kimm et al. (2017) Kimm T., Katz H., Haehnelt M., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS, 466, stx052
  * Kimm et al. (2018) Kimm T., Haehnelt M., Blaizot J., Katz H., Michel-Dansac L., Garel T., Rosdahl J., Teyssier R., 2018, MNRAS, 475, 4617
  * Koudmani et al. (2022) Koudmani S., Sijacki D., Smith M. C., 2022, MNRAS, 516, 2112
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Krumholz & Tan (2007) Krumholz M. R., Tan J. C., 2007, ApJ, 654, 304
  * Kulsrud & Pearce (1969) Kulsrud R., Pearce W. P., 1969, ApJ, 156, 445
  * Kuwabara & Ko (2020) Kuwabara T., Ko C.-M., 2020, ApJ, 899, 72
  * Leitherer et al. (1999) Leitherer C., et al., 1999, ApJS, 123, 3
  * Leroy et al. (2023) Leroy A. K., et al., 2023, ApJL, 944, L10
  * Lopez-Rodriguez et al. (2021) Lopez-Rodriguez E., Guerra J. A., Asgari-Targhi M., Schmelz J. T., 2021, ApJ, 914, 24
  * Maiolino & Mannucci (2019) Maiolino R., Mannucci F., 2019, A&ARv, 27, 3
  * Mannucci et al. (2010) Mannucci F., Cresci G., Maiolino R., Marconi A., Gnerucci A., 2010, MNRAS, 408, 2115
  * Mao & Ostriker (2018) Mao S. A., Ostriker E. C., 2018, ApJ, 854, 89
  * Marasco et al. (2023) Marasco A., et al., 2023, A&A, 670, A92
  * Martin-Alvarez et al. (2018) Martin-Alvarez S., Devriendt J., Slyz A., Teyssier R., 2018, MNRAS, 479, 3343
  * Martin-Alvarez et al. (2020) Martin-Alvarez S., Slyz A., Devriendt J., Gómez-Guijarro C., 2020, MNRAS, 495, 4475
  * Martin-Alvarez et al. (2021) Martin-Alvarez S., Katz H., Sijacki D., Devriendt J., Slyz A., 2021, MNRAS, 504, 2517
  * Martin-Alvarez et al. (2022) Martin-Alvarez S., Sijacki D., Haehnelt M. G., Farcy M., Dubois Y., Belokurov V., Rosdahl J., Lopez-Rodriguez E., 2022, MNRAS, 000, 1
  * Mashchenko et al. (2006) Mashchenko S., Couchman H. M. P., Wadsley J., 2006, Nature, 442, 539
  * McQuinn et al. (2019) McQuinn K. B. W., van Zee L., Skillman E. D., 2019, ApJ, 886, 74
  * Miyoshi & Kusano (2005) Miyoshi T., Kusano K., 2005, J. Comput. Phys., 208, 315
  * Mocz et al. (2016) Mocz P., Pakmor R., Springel V., Vogelsberger M., Marinacci F., Hernquist L., 2016, MNRAS, 463, 477
  * Morlino & Caprioli (2012) Morlino G., Caprioli D., 2012, in AIP Conference Proceedings. pp 241–244, doi:10.1063/1.4772242
  * Moster et al. (2013) Moster B. P., Naab T., White S. D. M., Moster B. P., Naab T., White S. D. M., 2013, MNRAS, 428, 3121
  * Mulcahy et al. (2014) Mulcahy D. D., et al., 2014, A&A, 568, A74
  * Murray et al. (2005) Murray N., Quataert E., Thompson T. A., 2005, ApJ, 618, 569
  * Naab & Ostriker (2017) Naab T., Ostriker J. P., 2017, ARA&A, 55, 59
  * Naab et al. (2007) Naab T., Johansson P. H., Ostriker J. P., Efstathiou G., Naab T., Johansson P. H., Ostriker J. P., Efstathiou G., 2007, ApJ, 658, 710
  * Nuñez-Castiñeyra et al. (2022) Nuñez-Castiñeyra A., Grenier I. A., Bournaud F., Dubois Y., Kamal Youssef F. R., Hennebelle P., 2022, arXiv e-prints, p. arXiv:2205.08163
  * Olesen (1993) Olesen P., 1993, arXiv e-prints, pp hep–ph/9310275
  * Pakmor & Springel (2013) Pakmor R., Springel V., 2013, Monthly Notices of the Royal Astronomical Society, 432, 176
  * Pakmor et al. (2016) Pakmor R., Pfrommer C., Simpson C. M., Springel V., 2016, ApJ, 824, L30
  * Peschken et al. (2022) Peschken N., Hanasz M., Naab T., Wóltański D., Gawryszczak A., 2022, arXiv e-prints, 000, arXiv:2210.17328
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 473, 4077
  * Powell et al. (2011) Powell L. C., Slyz A., Devriendt J., 2011, MNRAS, 414, 3671
  * Power et al. (2003) Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 338, 14
  * Ptuskin et al. (1997) Ptuskin V. S., Voelk H. J., Zirakashvili V. N., Breitschwerdt D., Ptuskin V. S., Voelk H. J., Zirakashvili V. N., Breitschwerdt D., 1997, A&A, 321, 434
  * Pudritz et al. (1989) Pudritz R. E., Silk J., Pudritz R. E., Silk J., 1989, ApJ, 342, 650
  * Quataert et al. (2021a) Quataert E., Jiang Y.-F., Thompson T. A., 2021a, MNRAS, 510, 920
  * Quataert et al. (2021b) Quataert E., Thompson T. A., Jiang Y.-F., 2021b, MNRAS, 510, 1184
  * Rasera & Teyssier (2006) Rasera Y., Teyssier R., 2006, A&A, 445, 1
  * Rey et al. (2023) Rey M. P., Katz H. B., Cameron A. J., Devriendt J., Slyz A., 2023, arXiv e-prints, p. arXiv:2302.08521
  * Rhee et al. (2022) Rhee J., Elahi P., Yi S. K., 2022, ApJ, 927, 129
  * Rieder & Teyssier (2017) Rieder M., Teyssier R., 2017, MNRAS, 471, 2674
  * Rodríguez Montero et al. (2022) Rodríguez Montero F., Martin-Alvarez S., Sijacki D., Slyz A., Devriendt J., Dubois Y., 2022, MNRAS, 511, 1247
  * Rosdahl & Teyssier (2015) Rosdahl J., Teyssier R., 2015, MNRAS, 449, 4380
  * Rosdahl et al. (2018) Rosdahl J., et al., 2018, MNRAS, 479, 994
  * Rosen & Bregman (1995) Rosen A., Bregman J. N., 1995, ApJ, 440, 634
  * Ruszkowski et al. (2017) Ruszkowski M., Yang H.-Y. K., Zweibel E., 2017, ApJ, 834, 208
  * Salem et al. (2014) Salem M., Bryan G. L., Hummels C., 2014, ApJ, 797, L18
  * Schlickeiser et al. (2009) Schlickeiser R., Schlickeiser R. 2009, MNRAS, 398, 1483
  * Schmidt (1959) Schmidt M., 1959, ApJ, 129, 243
  * Shadmehri (2009) Shadmehri M., 2009, MNRAS, 397, 1521
  * Silk (2010) Silk J., 2010, Proc. Int. Astron. Union, 6, 273
  * Skilling (1971) Skilling J., 1971, ApJ, 170, 265
  * Skilling (1975) Skilling J., 1975, MNRAS, 172, 557
  * Smith et al. (2018) Smith M. C., Sijacki D., Shen S., 2018, MNRAS, 478, 302
  * Smith et al. (2019) Smith M. C., Sijacki D., Shen S., 2019, MNRAS, 485, 3317
  * Socrates et al. (2006) Socrates A., Davis S. W., Ramirez-Ruiz E., 2006, ApJ, 687, 202
  * Somerville & Davé (2015) Somerville R. S., Davé R., 2015, ARA&A, 53, 51
  * Springel & Hernquist (2003) Springel V., Hernquist L., 2003, MNRAS, 339, 312
  * Stinson et al. (2006) Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn T., 2006, MNRAS, 373, 1074
  * 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
  * Summers et al. (2003) Summers L. K., Stevens I. R., Strickland D. K., Heckman T. M., 2003, MNRAS, 342, 690
  * Tasker (2011) Tasker E. J., 2011, ApJ, 730, 11
  * Teyssier (2002) Teyssier R., 2002, A&A, 385, 337
  * Teyssier et al. (2006) Teyssier R., Fromang S., Dormy E., 2006, J. Comput. Phys., 218, 44
  * Teyssier et al. (2013) Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MNRAS, 429, 3068
  * Thilker et al. (2023) Thilker D. A., et al., 2023, ApJL, 944, L13
  * Thornton et al. (1998) Thornton K., Gaudlitz M., Janka H., Steinmetz M., 1998, ApJ, 500, 95
  * Tillson et al. (2015) Tillson H., Devriendt J., Slyz A., Miller L., Pichon C., 2015, MNRAS, 449, 4363
  * Tortora et al. (2022) Tortora C., Hunt L. K., Ginolfi M., 2022, A&A, 657, A19
  * Tóth et al. (2000) Tóth G., Tóth Gábor 2000, JCoPh, 161, 605
  * Tremonti et al. (2004) Tremonti C. A., et al., 2004, ApJ, 613, 898
  * Tweed et al. (2009) Tweed D., Devriendt J., Blaizot J., Colombi S., Slyz A., 2009, A&A, 506, 647
  * Vachaspati (1991) Vachaspati T., 1991, Phys. Lett. B, 265, 258
  * Vachaspati (2021) Vachaspati T., 2021, Rep. Prog. Phys., 84, 074901
  * Vazza et al. (2017) Vazza F., Jones T. W., Brüggen M., Brunetti G., Gheller C., Porter D., Ryu D., 2017, MNRAS, 464, 210
  * Veilleux et al. (2005) Veilleux S., Cecil G., Bland-Hawthorn J., 2005, ARA&A, 43, 769
  * Vogelsberger et al. (2020) Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nat. Rev. Phys., 2, 42
  * Wadepuhl & Springel (2011) Wadepuhl M., Springel V., 2011, MNRAS, 410, 1975
  * Wagner et al. (2005) Wagner A. Y., Falle S. A. E. G., Hartquist T. W., Pittard J. M., 2005, A&A, 430, 567
  * Werhahn et al. (2021) Werhahn M., Pfrommer C., Girichidis P., 2021, MNRAS, 508, 4072
  * Whalen et al. (2008) Whalen D., van Veelen B., O’Shea B. W., Norman M. L., 2008, ApJ, 682, 49
  * Williams et al. (1997) Williams J. P., McKee C. F., Williams J. P., McKee C. F., 1997, ApJ, 476, 166
  * Wise et al. (2012) Wise J. H., Turk M. J., Norman M. L., Abel T., 2012, ApJ, 745, 50
  * Wolfire et al. (1995) Wolfire M. G., et al., 1995, ApJ, 443, 152
  * Xu & Lazarian (2022) Xu S., Lazarian A., 2022, ApJ, 927, 94
  * Zubko et al. (2004) Zubko V., Dwek E., Arendt R. G., 2004, ApJS, 152, 211
  * Zuckerman & Evans (1974) Zuckerman B., Evans N. J. I., 1974, ApJ, 192, L149
  * Zweibel (2017) Zweibel E. G., 2017, Phys. Plasmas, 24, 055402
  * de Blok (2010) de Blok W. J. G., 2010, Adv. Astron., 2010, 1