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

Francisco Rodríguez Montero, Sergio Martin-Alvarez, Adrianne Slyz, Julien

Devriendt, Yohan Dubois, Debora Sijacki

► To cite this version:

Francisco Rodríguez Montero, Sergio Martin-Alvarez, Adrianne Slyz, Julien Devriendt, Yohan Dubois, et al.. The impact of cosmic rays on the interstellar medium and galactic outflows of Milky Way analogues. Monthly Notices of the Royal Astronomical Society, 2024, 530 (4), pp.3617-3640. 10.1093/mn-ras/stae1083. hal-04179299

HAL Id: hal-04179299

https://hal.science/hal-04179299v1

Submitted on 26 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The impact of cosmic rays on the interstellar medium and galactic outflows of Milky Way analogues

Francisco Rodríguez Montero<sup>1,★</sup> Sergio Martin-Alvarez<sup>1,2</sup> Adrianne Slyz,<sup>1</sup> Julien Devriendt,<sup>1,3</sup> Yohan Dubois<sup>4</sup> and Debora Sijacki<sup>5</sup>

Accepted 2024 April 16. Received 2024 March 22; in original form 2023 July 25

## 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, while 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 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 ~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 supernova (SN) events exploding at lower densities, where they have a higher impact. SN-injected CRs produce enhanced, multiphase 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 ionized gas, warm neutral and cold gas phases contribute significantly at high redshifts. Importantly, our work shows that future *JWST* observations of galaxies and their multiphase outflows across cosmic time have the ability to constrain the role of CRs in regulating star formation.

**Key words:** methods: numerical – stars: formation – cosmic rays – ISM: jets and outflows – galaxies: evolution – galaxies: formation – galaxies: spiral.

## 1 INTRODUCTION

Realistically, modelling galaxy formation within the Lambda cold dark matter cosmological paradigm has proven to be an extremely complex task, as its multiscale and multiphysics 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; Silk 2010; Dobbs, Burkert & Pringle 2011; Hopkins, Quataert & Murray 2011; Tasker 2011; 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. Springel & Hernquist 2003; Stinson et al. 2006; Teyssier et al. 2013), including such a considerable energy injection through SN or active galactic nuclei (AGNs) 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, Hunt & Ginolfi 2022), the metal enrichment of the circum- (CGM) and intergalactic medium (IGM; e.g. Veilleux, Cecil & Bland-Hawthorn 2005; Brooks et al. 2007; Chisholm, Tremonti & Leitherer 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, 2016), the flattening of inner dark matter (DM) profiles (e.g. Mashchenko, Couchman & Wadsley 2006; de Blok 2010; Governato et al. 2010; Teyssier et al. 2013) and the formation of bulge-less

<footnote n="1"><sup>1</sup>Sub-department of Astrophysics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="2"><sup>2</sup>Kavli Institute for Particle Astrophysics and Cosmology (KIPAC), Stanford University, Stanford, CA 94305, USA</footnote>

<footnote n="3"><sup>3</sup>CNRS, Centre de Recherche Astrophysique de Lyon, Université de Lyon, Université Lyon 1, ENS de Lyon, UMR 5574, F-69230 Saint-Genis-Laval, France</footnote>

<footnote n="4"><sup>4</sup>Institut d'Astrophysique de Paris, UMR 7095, Sorbonne Université, CNRS, 98 bis boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="5"><sup>5</sup>Institute of Astronomy and Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote>*E-mail: francisco.rodriguezmontero@physics.ox.ac.uk</footnote>
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 multiphase structure of the ISM is beginning to be spatially resolved. This has driven the development of more physically motivated treatments of SN 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 photoionization heating and radiation pressure (Murray, Quataert & Thompson 2005; Hopkins et al. 2011; Agertz et al. 2013; Rosdahl & Teyssier 2015; Emerick, Bryan & Low 2018; Smith, Sijacki & Shen 2018, 2019), stellar winds (Hopkins et al. 2011; Agertz et al. 2013; Geen et al. 2015; Fierlinger et al. 2016; Smith, Sijacki & Shen 2019), or direct pressure from Ly$\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 remain 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 (~$10^6$ K) gas (Heckman et al. 1995; Summers et al. 2003; McQuinn, van Zee & Skillman 2019; Concas et al. 2022; Marasco et al. 2023), 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 that are accelerated by dif-
fusive shock acceleration (e.g. Axford 1981; Bell 1978; Blandford &
Eichler 1987), act as a non-thermal source of energy, and pioneering
numerical investigations (e.g. Jubelgas et al. 2008; Wadepuhl &
Springel 2011; Booth et al. 2013; Salem, Bryan & Hummels 2014)
of their effects indicate that they can have a significant impact on star
formation. Observations of the Milky Way (MW) and nearby galaxies
suggest that CRs and magnetic fields are in equipartition with the
thermal and turbulent energies in the disc mid-plane (Boulares &
Cox 1990; Basu & Roy 2013; Beck 2016; Zweibel 2017; Lopez-
Rodriguez et al. 2021). Furthermore, the peak of their spectrum at
~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, McKenzie & Voelk 1991; Ptuskin
et al. 1997; Socrates, Davis & Ramirez-Ruiz 2006; Everett et al.
2008; Mao & Ostriker 2018; Quataert, Jiang & Thompson 2021a;
Quataert, Thompson & Jiang 2021b). These ~1 GeV CRs behave as
a relativistic fluid with a soft equation of state (i.e. adiabatic index
$\gamma_{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 (~$10^{6-7}$ Myr)
via Coulomb and hadronic interactions (Ensslin 2004; Schlickeiser

& Schlickeiser 2009) than the radiative cooling of thermal gas. They rapidly diffuse away from their injection sites as their diffusion lengths are larger than typical molecular clouds (~100 pc), reaching the more dilute gas immediately above the galactic disc, 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, Strong & Girichidis 2021, and references therein): CRs thicken the gaseous disc by altering the pressure-gravity balance, changing the morphology of the galaxy (Booth et al. 2013; Salem et al. 2014; Ruszkowski, Yang & Zweibel 2017; Buck et al. 2020; Dashyan & Dubois 2020; Chan et al. 2022; Farcy et al. 2022; Nuñez-Castiñeyra et al. 2022; Martin-Alvarez et al. 2022). In addition, CR feedback significantly affects the star formation history (SFH) of galaxies by 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 launching more heavily mass-loaded, denser and colder outflows (e.g. Booth et al. 2013; Hanasz et al. 2013; Hopkins et al. 2021a; Farcy et al. 2022; Peschken et al. 2023).

This modelling progress has established that CRs in the ~1 GeV regime likely constitute an important ingredient of galaxy formation. However, a clear understanding of how they couple with other physical processes of relevance remains far from settled. 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, which can only be captured by fully self-consistent cosmological simulations (see e.g. Ji et al. 2020, 2021; Beckmann et al. 2022; Martin-Alvarez et al. 2022).

Although ~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 magnetohydrodynamic (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, Ostriker & Jiang 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 magnetothermal 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, Marcowith & Dubois 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 multiphase outflows (Salem et al. 2014; Girichidis et al. 2016; Kim & Ostriker 2018; Hopkins et al. 2021a; Martin-Alvarez et al. 2022) which require to be carefully studied in an explicit cosmological context.
This is the first paper in a series in which we study the effects of CR feedback in galaxies with a resolved multiphase 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 MW-like galaxy (called NUT) formation simulation set-up, and our CR and MHD modelling. We then discuss the global effects of magnetic fields and CRs in the evolution of NUT in Section 3.1, demonstrating how 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 summarize our main results, discuss their limitations and the prospect for future studies of CR feedback.

# 2 NUMERICAL METHODS

## 2.1 The RAMSES code

The full suite of new simulations presented in this work has been generated using the RAMSES code<fnref n="1" /> (Teyssier 2002), an Adaptive Mesh Refinement (AMR) Eulerian code with an octree grid structure. RAMSES uses collisionless particles to model the evolution of the DM and stellar components. These are coupled to each other and to the gas component through a gravitational potential that obeys the Poisson equation and is solved on the adaptive mesh. In order to evolve in time the MHD properties 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 (Fromang, Hennebelle & Teyssier 2006; Teyssier, Fromang & Dormy 2006). This method ensures that the solenoidal condition (i.e. $\vec{\nabla} \cdot \vec{B} = 0$) is fulfilled to numerical machine precision, thus preventing numerical spurious growth of magnetic field and cell noise (see Tóth, Tóth & Gábor 2000; Mocz et al. 2016, for a comparison of CT with other MHD methods).

### 2.1.1 Cosmic ray MHD

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:

$$ \frac{\partial \rho}{\partial t} + \vec{\nabla} \cdot (\rho \vec{u}) = 0, \qquad (1) $$

$$ \frac{\partial \rho \vec{u}}{\partial t} + \vec{\nabla} \cdot \left( \rho \vec{u} \vec{u} + P_{\text{tot}} - \frac{\vec{B} \cdot \vec{B}}{4\pi} \right) = \rho \vec{g}, \qquad (2) $$

$$ \frac{\partial e}{\partial t} + \vec{\nabla} \cdot \left( (e+P)\vec{u} - \frac{\vec{B}(\vec{B} \cdot \vec{u})}{4\pi} \right) = \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}}, \qquad (3) $$

$$ \frac{\partial \vec{B}}{\partial t} - \vec{\nabla} \times (\vec{u} \times \vec{B}) = 0, \qquad (4) $$

$$ \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}}) = -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}}, \qquad (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_{\text{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_{\text{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_{\text{CR}} = 4/3$ for an ideal non-relativistic and relativistic gas, respectively. This modelling approximation can be adopted for the evolution of low-energy (~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}_{\text{CR,d}} = -\kappa \vec{b}(\vec{b} \cdot \vec{\nabla} e_{\text{CR}})$ the anisotropic diffusion flux (κ is the diffusion coefficient, assumed constant in this work), $\mathcal{L}_{\text{rad}} = \mathcal{L}_{\text{rad,th}} + \mathcal{L}_{\text{rad,CR} \to \text{th}}$ the total radiative losses which includes thermal and CR radiative losses, respectively.<fnref n="2" /> The CR loss term encompasses both energy losses due to hadronic and Coulomb interactions $\mathcal{L}_{\text{CR}} = -7.51 \times 10^{-16}(n_e/\text{cm}^{-3})(e_{\text{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}_{\text{CR} \to \text{th}} = 2.63 \times 10^{-16}(n_e/\text{cm}^{-3})(e_{\text{CR}}/\text{erg cm}^{-3})\text{erg s}^{-1}$.

In the CR energy equation (equation 5), $\vec{F}_{\text{st}} = (e_{\text{CR}} + P_{\text{CR}})\vec{u}_{\text{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; Thomas, Pfrommer & Pakmor 2023), allowing for streaming velocities exceeding $\vec{u}_{\text{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; Thomas et al. 2023) for future work. Furthermore, since CRs scatter off the Alfvén waves, they experience a drag force, transferring energy to the thermal gas at a rate $\mathcal{L}_{\text{st}} = -\text{sign}(\vec{b} \cdot \vec{\nabla} e_{\text{CR}})\vec{u}_{\text{A}} \cdot \vec{\nabla} P_{\text{CR}}$. The typical value of $\kappa \simeq 3 \times 10^{28}(R_{\text{CR}}/\text{GV})^{0.34} \text{cm}^2 \text{s}^{-1}$ depends on the particle rigidity $R_{\text{CR}}$ (Jóhannesson, Porter & Moskalenko 2019), which for 1 GeV proton has a value of $R_{\text{CR}} = 1$ GV. Previous efforts to constrain $\kappa$ in CR simulations by comparing Fermi LAT γ-ray observations (e.g. FIRE-2, AREPO, and RAMSES; Chan et al. 2019; Werhahn, Pfrommer & Girichidis 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 γ-ray emissivity per atom can reduce the discrepancy by factors of up to 2.5, these are not sufficient to bring the different simulations in agreement. They suggest that discrepancies in the thermal state of the gas and magnetic field properties, when evolved with different MHD codes, more likely source these differences. The exploration of how different CR transport mechanisms and diffusion coefficients affect the γ-ray luminosity is the scope of an upcoming paper (Rodríguez Montero et al., in preparation). We adopt here a value of $\kappa = 3 \times 10^{28} \text{cm}^2 \text{s}^{-1}$ favoured by the results of isolated CR simulations using versions of RAMSES similar to the one used in this work (Dashyan & Dubois 2020; Farcy et al. 2022; Nuñez-Castiñeyra et al. 2022), and leave the

<footnote n="1">¹https://bitbucket.org/rteyssie/ramses/src/master/</footnote>

<footnote n="2">²The modelling of radiative losses is described in Section 2.2.</footnote>
<table><caption>**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, the employed solver and the details of the CR physics if CRs are present. All simulations use the same initial conditions, an MTT star formation prescription and the canonical configuration for mechanical SN feedback. 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.</caption><thead><tr><th>Simulation</th><th>Solver</th><th>CR physics</th></tr></thead><tbody><tr><td>HD</td><td>Hydro</td><td><i>x</i></td></tr><tr><td>MHD</td><td>MHD</td><td><i>x</i></td></tr><tr><td>CRMHD</td><td>CRMHD</td><td>Streaming and <i>&kappa;</i> = 3 &times; 10<sup>28</sup> cm<sup>2</sup> s<sup>-1</sup></td></tr></tbody></table>

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, Slyz & Devriendt 2011; Geen, Slyz & Devriendt 2013) generated at $z = 499$ with cosmological parameters consistent with the results of WMAP5 (matter density $\Omega_m = 0.258$, dark energy density $\Omega_\Lambda = 0.742$, baryon density $\Omega_b = 0.045$, amplitude of the linear

power spectrum at the 8 Mpc scale $\sigma_8 = 0.8$ and Hubble constant $H_0 = 72 \text{ km s}^{-1} \text{ 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 the high resolution calculation is performed. Within this zoom-in region a $M_{\text{vir}}(z=0) \simeq 5 \times 10^{11} \text{ M}_\odot$ DM halo hosts an MW-like galaxy (from this point simply identified as the NUT galaxy) with a DM particle resolution of $m_{\text{DM}} = 5.5 \times 10^4 \text{ M}_\odot$. To probe the multiphase 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 DM) 8 times larger than $m_{\text{DM}}\Omega_b/\Omega_m$ are subdivided into 8 smaller cells of equal size, down to a minimum cell size $\Delta x_{\min} \simeq 23 \text{ pc}$ approximately constant in proper units (an extra level of refinement is triggered for expansion factors of $a_{\text{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 lives 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 appear 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

<figure><img src="image_2.png" /><figcaption>**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 a 600 kpc projection 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 disc 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 outflow.</figcaption></figure>
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, Dwek & Arendt (2004). At $z = 1.5$, NUT shows signatures of a large feedback event: hot bipolar outflows emerge perpendicular to the cold, dense and magnetized disc. These outflows clear denser and colder gas present in the CGM, creating bubbles seen in the gas density, and both magnetic and CR energies. The face-on view through JWST filters shows a late-type spiral with a nuclear cluster at its centre dominated by an old stellar population, and a disc 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 ~$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_H = 0.01 \text{ cm}^{-3}$ considered to be self-shielded from UV radiation with an exponential dependence on $n_H$. Furthermore, despite the high resolution of our simulation, we lack the cooling modelling capabilities to resolve primordial (Population III) SN explosions in minihaloes (e.g. Whalen et al. 2008; Wise et al. 2012) that would potentially raise the initial gas metallicity to ~$10^{-3} Z_\odot$ and hence increase cooling at high redshift. To avoid this issue, we initialize 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 magnetohydro-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 \text{ cm}^{-3}$ and the gravitational pull is larger than the combined support from turbulent, thermal, magnetic, and CR pressure.<fnref n="3" /> We determine the rate of gas to stellar mass conversion following a Schmidt law (Schmidt 1959):

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

where we allow $\epsilon_{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) [generalized 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 simply 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 an 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 an efficiency of $\epsilon_{SN} = 10^{51}$ erg/(10 $M_\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 chromodynamics phase transitions Vachaspati 1991; Hollenstein et al. 2008) 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 ~$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-axis.<fnref n="4" /> 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 initialization directly provides the expected strengths of $1–10 \mu\text{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_{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 OZYMANDIAS<fnref n="5" />, a new cosmological simulation analysis package that 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 parallelized routines to access the raw RAMSES data. It takes a cosmological simulation snap-

<footnote n="3">³We have monitored 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 ~10–30 per cent, although in rare occasions the support of CR pressure can be enough to lower this efficiency to 0 per cent.</footnote>

<footnote n="4">⁴In 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 magnetization of the large scale environment changes.</footnote>

<footnote n="5">⁵https://github.com/Currodri/ozymandias</footnote>
shot and outputs a convenient galaxy and halo catalogue following the 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 per cent 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 three-dimensional (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, Elahi & Yi (2022), a robust, scale-free, phase-space structure-finder. Once haloes and galaxies are identified, a hierarchical structure is built between them, linking galaxies with their host haloes 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 at any previous snapshot. We use this package to identify the main galaxy, and refer to the galactic region as the sphere with radius $0.2R_{vir,DM}$, $R_{vir,DM}$ being the host DM halo virial radius. Quantities such as the galaxy stellar mass $M_*$ and gas mass $M_{gas}$ are measured within this spherical volume at every snapshot.

# 3 RESULTS

## 3.1 The global impact of cosmic rays on the evolution of NUT

In Fig. 2, we show density-weighted projections of NUT at $z = 1.5$ for the three simulations [hydrodynamical (HD), MHD, and cosmic ray MHD (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 disc 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 disc galaxy within a virialized halo: a cold ($T < 10^4$ K) thin disc with close to solar metallicity $Z_\odot$ 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 disc can still be clearly seen in the projections, there are significant differences with respect to the no-CR simulations. First, the gas disc 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 (resulting in a CR pressure-dominated halo as seen in Appendix A). Against this lower temperature ‘background’, recent hot outflows can clearly

be distinguished in the edge-on view, although the lower average metallicity of the disc 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_\odot$.

Furthermore, the mock JWST images of the HD and MHD simulations show NUT as a disc 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 disc instabilities persist due to inefficient feedback (e.g. Bournaud et al. 2013; Inoue et al. 2016). 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 disc, small star-forming regions are seen all across the disc, 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):

$$ \log \Sigma_{1.5} = \log(M_*/M_\odot) - 1.5 \log(r_{50}/\text{kpc}), \quad (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.

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-hand panel) and gas (right-hand 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} M_\odot$ at $\sim 1.5$ 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} M_\odot$ at $z = 1.5$, in agreement with the equivalent mechanical feedback with real delay (MFBm) model in Kimm et al. (2015). The evolution of $M_{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 time-scale 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 (e.g. Martin-Alvarez et al. 2018). Fig. 3 includes the MHD run with a primordial magnetic field of $3 \times 10^{-12} \text{ G}$ but no cosmic rays. This illustrates that the moderate primordial magnetization employed here does not significantly affect the global properties of the galaxy (Su et al. 2017; Pillepich et al. 2018; Martin-Alvarez et al. 2018; Hopkins et al. 2020).

However, the addition of CRs gives rise to a radical difference in the SFH of the NUT galaxy: the steepness of the stellar mass growth is reduced during the accretion phase, resulting in $M_* \sim 2 \times 10^9 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 only exhibits a $\sim 2-4$ deficit in $M_{gas}$ inside the galactic region. Therefore, the decrease in $M_*$ cannot be

<footnote n="6">⁶https://www.hdfgroup.org/solutions/hdf5/</footnote>
<figure><img src="image_3.png" /><figcaption>**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_{vir,DM}$ and galaxy compactness $\log \Sigma_{1.5}$ (Barro et al. 2013) are listed, with a larger value indicating a more compact stellar emission. Compared to HD and MHD simulations, the disc 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 disc. In this case, star-forming regions are scattered across spiral arms and clumps all over the disc.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Evolution of the stellar (left-hand panel) and gas (right-hand panel) masses contained within the galactic region (i.e. $<0.2R_{vir,DM}$) of the NUT galaxy in the HD (blue line), MHD (purple line), and CRMHD (green line) simulations. The right-hand 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, the stellar mass remains 4 times lower.</figcaption></figure>

solely attributed to the depletion of $M_{\text{gas}}$ from the galactic region. We also include on the right-hand panel the cold gas mass (dashed lines) as defined by the constant entropy limit<fnref n="7" /> 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$. These 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 undergoes 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_{\text{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 gas and DM mass enclosed within $R_{vir,DM}$ plus 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, Wechsler & Conroy (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 upper limit of total conversion of gas into stars for a given halo, given by $f_b M_{\text{halo}}$, where $f_b$ is the universal baryon fraction defined as $\Omega_b/\Omega_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 DM halo mass. The HD run begins with $\sim 32$ per cent baryon conversion efficiency, which grows to $\sim 56$ per cent by $z = 2$, while the MHD run reaches a slightly lower efficiency of $\sim 47$ per cent. The CRMHD run shows a consistent decrease of the baryonic conversion efficiency across all redshifts, with a maximum value of $\sim 14$ per cent at $z = 2$ and a minimum of 5 per cent 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 small haloes with masses of $10^{10} M_{\odot}$ or lower between $z = 6$ and $z = 8$ are scarce, increasing the uncertainty of abundance matching models at high redshift. Additionally, the observational stellar mass functions used in these works are subject to observational incompleteness below $10^7 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.

<footnote n="7">⁷The cold phase specific entropy limit is given by $s \le 4.4 \times 10^8 \text{ erg K}^{-1} \text{ g}^{-1}$, e.g. corresponding to a density of $7.8 \times 10^{-23} \text{ g cm}^{-3}$ at 300 K.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 4. Stellar mass to halo mass function of the NUT galaxy for our fiducial models: HD (blue), MHD (purple), 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 to transform the majority of their baryonic content into stars, while the CRMHD run from early on ($z \sim 8$) sits ~1 dex below, approaching the prediction by Behroozi et al. (2013) at $z \sim 1.5$.</figcaption></figure>

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

As shown in Section 3.1, the CRMHD simulation contains a factor of ~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 dynamical time of the galaxy at $z \sim 2.5$, weighted by their distance in time to $z \sim 2.5$. We estimate the dynamical time of the galaxy $\tau_{dyn}$ as the time required for a test particle to complete one full orbit at $r = 0.2R_{vir, DM}$ with circular velocity:

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

where $M(<r) = M_b + M_{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_{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-hand panel of Fig. 5 shows that the low density gas ($10^{-26} \text{ g cm}^{-3}$ to $10^{24} \text{ g cm}^{-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

<figure><img src="image_6.png" /><figcaption>Figure 5. Stacked gas temperature-density phase diagrams of the NUT galaxy ($<0.2R_{vir, DM}$) in HD (left-hand panel), MHD (middle panel) and CRMHD (right-hand 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 enclose regions where ~80 per cent of the gas with $\epsilon_{SF}^{\text{MTT}} \ge 1$ per cent lie. The less intense white hue within the black contour of the CRMHD run indicates a much lower total gas mass in the regime favourable to star formation.</figcaption></figure>
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), Martin-Alvarez et al. (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 retaining their temperature. This is reflected in the size of the 'dark' regions of the phase diagram with densities below $10^{-24}$ 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 per cent in HD to 55 per cent in MHD). The hot phase only contributes 3 per cent to the total gas mass in both HD and MHD.

However, the gas phase diagram for the CRMHD run (the right-hand panel of Fig. 5) shows marked differences compared with the non-CR models: (a) at densities $\rho < 10^{-24}$ g cm$^{-3}$, gas is no longer confined to temperatures above $10^4$ K when $\rho < 10^{-24}$ 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}$ 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 2 compared to the HD model as this gas now resides in the warm medium (as indicated by the less intense white hue visible in this region on the right-hand panel of Fig. 5). Notably, it is in these dense regions of the phase diagram that we expect gas to undergo a more efficient star formation. Overlaid black contours enclose the regions in temperature-density space that contain 80 per cent of the star-forming gas with efficiency (as given by the MTT model) above 1 per cent. It is readily discernible that the gradual incorporation of additional sources of support (such as magnetic and/or CR pressures) results in a reduction of the area encompassed by these contours. While the HD and MHD star-forming gas contributes to 3 per cent and 2 per cent of the gas in the galactic region, for CRMHD is limited to 0.06 per cent. This considerable reduction of star-forming 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 disc 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 $10^{-24}$ g cm$^{-3}$ to $10^{-22}$ g cm$^{-3}$ in the CRMHD run creates a more pronounced concave inflection point towards lower $T$ at $\rho \sim 10^{-23}$ g cm$^{-3}$; (ii) gas mass

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Evolution of specific thermal (solid), magnetic (dotted), and CR (dashed) energies within the galactic region ($0.2 R_{\text{vir,DM}}$) of NUT. We compare HD (blue line), MHD (purple lines), and CRMHD (green lines) runs. Thermal specific energy shows a decrease of a factor of ~3–4 between high and low redshift, with the CRMHD specific thermal energy consistently higher than that of the other two runs down to $z \sim 2$. From very early on, the CR specific energy clearly dominates (~1–2 dex above thermal specific energy) the energy budget in the CRMHD simulation. Magnetic specific energy plays a secondary role, being well below thermal or CR specific energies at all times and in all runs. Both CRMHD and MHD run exhibit a weak increase in magnetic specific energy with cosmic time.</figcaption></figure>

accumulates more tightly around the thermal instability shoulder at densities below $10^{-24}$ g cm$^{-3}$, driving the median line closer to this feature of the phase diagram as the relative amount of gas above $T \ge 10^5$ K is greatly reduced, particularly at higher densities.

To illustrate the importance of CRs in energetic terms, Fig. 6 shows the evolution of the density-weighted average specific energies within the $0.2 R_{\text{vir,DM}}$ of NUT for the HD, MHD, and CRMHD simulations. This provides us with an intuition about which energy components will be most influential in driving the evolution of the denser (star-forming) regions of the galaxy. The thermal specific energy only shows a mild decrease with cosmic time in all three runs, with the thermal energy somewhat larger in the CRMHD run than in the HD and MHD simulations down to $z \sim 2$. This larger density-weighted average is likely related to the warm-diffuse phase being more mass dominant in the CRMHD simulation. Both the MHD and CRMHD runs display sub-dominant magnetic specific energies with a minor growth trend (more pronounced in the CRMHD simulation) over the entire redshift range. Although this may indicate extra magnetic field amplification taking place within the galactic region, a careful examination of magnetic field properties and amplification is beyond the scope of the current study. Finally, when present, the CR specific energy clearly dominates the energy budget, exceeding the thermal specific energy by up to ~2 dex at times. This supports the view that CRs play a dominant dynamical role, particularly in the dense gas within the galactic region. It is only at the very high densities of the galactic centre where thermal specific energy becomes comparable or even dominates over CR specific energy (see Appendix A).

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

### 3.3.1 Distribution of SN sites and feedback effectiveness

The modifications of the gas distribution across ISM phases also impact the distribution of SN explosions. In turn, these local
<figure><img src="image_8.png" /><figcaption>**Figure 7.** 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 ($10^{-22}$ g cm$^{-3}$ to $10^{-21}$ g cm$^{-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 10^{-27}$ g cm$^{-3}$ and a peak at $\sim 10^{-23}$ g cm$^{-3}$. SN events taking place at lower densities during this period are more effective in affecting the surrounding gas. This bimodality 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 becoming more resembling to that measured in the HD and MHD simulations.</figcaption></figure>

properties affect the efficiency of SN feedback in self-regulating star formation and modifying the ISM properties. We explore this issue by recording local gas properties in the central cell where individual SN explosions occur. In Fig. 7, we plot the distribution of gas densities at SN locations 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, characterizing 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 SN<fnref n="8" /> 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 10^{-27}$ g cm$^{-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 superbubbles.

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 (see the green dashed lines in the right-hand panel of Fig. 3), with the high-density peak of the bimodal distribution now at $\sim 10^{-23}$ g cm$^{-3}$ in the CRMHD run instead of $\sim 10^{-21}$ g cm$^{-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 ($\le 10^{-22}$ g cm$^{-3}$), with a large number of SNe exploding at $\rho \sim 10^{-27}$ g cm$^{-3}$. We attribute this to star formation happening in a less concentrated manner, throughout the galactic region, with subsequent 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. 7 shows gas density probability density functions (PDFs) for SN sites at lower redshifts (between $z=3$ and $z=2$). In this lower redshift range, we find that the CRMHD PDF begins to resemble to the ones obtained for the HD and MHD runs, although maintaining a higher fraction of SNe at lower densities (i.e. $\sim 10^{-27}$ g cm$^{-3}$). We know that NUT transitions to a rotationally supported (i.e. ratio of rotational velocity to dispersion velocity larger than unity) disc at $z \sim 3.5$ (Martin-Alvarez et al. 2020). We thus attribute the change in PDF bimodality primarily to a change in the geometrical morphology of the galaxy in the CRMHD simulation. At high redshift, the galaxy exhibits an irregular, near spherical gaseous density distribution. This means that cosmic rays emitted by an SN will have to travel a significant distance through dense, cold gas in all directions before they can escape. On the other hand, when the galactic disc becomes established at lower redshift, CR emanating from a SN going off within the thin disc will rapidly encounter diffuse gas in directions quasi-perpendicular to the disc. In other words, in a disc geometry the larger surface-to-mass ratio leads generally to CRs interacting less with the star-forming gas before they reach the edge of the galaxy. To provide an intuition of whether CR transport along

<footnote n="8">⁸At high densities the Sedov–Taylor phase will hardly be captured by our ~20 pc resolution, so the mechanical feedback method from Kimm et al. (2017) will mainly inject the momentum reached at the end of the snowplough phase.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Density-weighted magnetic field strength for the CRMHD run along the three cylindrical coordinates $B_r$ (disc radial; dotted red line), $B_\phi$ (toroidal; dashed dark-gold line), and $B_z$ (disc altitude; dot-dashed olive line), as well as the total field strength $|B|$.</figcaption></figure>

magnetic field lines facilitate such scenario, in Fig. 8 we compare the density-weighted average magnitude of the total magnetic field strength (green solid line) to its cylindrical components $B_z$ (dot-dashed olive line), $B_\phi$ (dashed dark golden line), and $B_r$ (red dotted line). Up until $z \sim 3.5$, all components have similar magnitudes, due to the spherical-like geometry of the system. During this period, CR transport along magnetic field lines is effectively isotropic. Below this redshift, the galactic disc formation is reflected in the growth of the toroidal component $B_\phi$ which becomes significantly larger than $B_z$ and $B_r$. Later on, the disturbance caused by a merger with a massive companion at $z \sim 2.3$ brings the three cylindrical magnetic components in equipartition once more before $B_\phi$ starts growing again. A dominant $B_\phi$ once the galaxy has established a rotationally supported disc implies that CR will preferentially follow magnetic field lines within the plane of the galaxy. However, the presence of a non-negligible $B_z$ facilitates the transport of CRs to higher disc altitudes, and due to the thin scale height of the disc compared with its radial dimension (scale-height $\sim 0.5$ kpc at $z \sim 2$, see also Martin-Alvarez et al. 2020), lead to their subsequent escape in the direction perpendicular to the disc. As this is the case here, and the $B_z$ component of the magnetic field remains significant compared to $B_\phi$ during the period when the galaxy has a disc morphology ($z \lesssim 3.5$), a large fraction of CRs will escape perpendicularly to the galactic disc plane.

Finally, although we do not aim to investigate in detail the properties and evolution of the magnetic field, the total magnetic field strength shown in Fig. 8 exhibits a certain amount of growth over the whole redshift range. This may be related to the CR-driven growth of magnetic fields reported by other works (Hanasz, Wóltański & Kowalik 2009; Siejkowski et al. 2009; Kulpa-Dybeł et al. 2015; Nuñez-Castiñeyra et al. 2022).

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

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

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

where we sum over all cells with positive radial velocities $v_{r,i}$ in a thin spherical shell of width $\delta r_{shell}$ at a distance $r$ from the galactic centre. We employ a shell thickness of 0.02 $R_{vir,DM}$, which provides a consistent measurement across simulations and cosmic time as it naturally accounts for galaxy growth (see Appendix B for a comparison of different shell widths). Each cell is characterized 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$ per cent) 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 outwards 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-haloes containing more than 1000 particles. The tidal radius is approximated as the Jacobi radius of a DM halo with the density profile of an isothermal sphere (see equation 8.107 of Binney & Tremaine 2008). Finally, we divide outflowing gas into separate phases as follows:

(i) Hot gas: $T > 10^5$ K,

(ii) Warm ionized gas: $10^5 > T > 9000$ K,

(iii) Warm neutral gas: $9000 > T > 1000$ K and

(iv) Cold gas: $T < 1000$ K.

In Fig. 9, we show the results of this outflow analysis when the thin spherical shell is located at a radial distance of 0.2 $R_{vir,DM}$ (evolving from $\sim 1$ kpc at $z = 10$ to $\sim 20$ kpc at $z = 1.5$). From top to bottom, we show $\dot{M}_{outflow}$ (top panel), star formation rate (SFR) averaged over 100 Myr (middle panel), and the mass loading factor $\eta = \dot{M}_{outflow}/\text{SFR}$ (bottom panel). Solid lines denote the running median calculated over 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_{vir,DM}$. We distinguish three distinct epochs when comparing the HD, MHD, and CRMHD simulations: (i) between $z \sim 13$ and $z \sim 7$ outflow rates in the CRMHD run reach $2 - 4$ M$_\odot$ yr$^{-1}$, whereas the HD and MHD runs are closer to $0.1 - 0.4$ M$_\odot$ yr$^{-1}$; (ii) between $z \sim 7$ and $z \sim 3$ outflow rates in the HD and MHD runs hover around $1$ M$_\odot$ yr$^{-1}$ (with the MHD outflows having overall larger values than their HD counterparts), corresponding to strong starburst episodes (however, given average SFRs of $10$ M$_\odot$ yr$^{-1}$, mass loading factors remain well below unity, while they consistently exceed $\eta \sim 1$ in the CRMHD run); and (iii) between $z \sim 3$ and $z \sim 1.5$, when the SFR in all three runs peaks around $10$ M$_\odot$ yr$^{-1}$ at $z \sim 2.5$ (and subsequently slowly declines to $2 - 4$ M$_\odot$ yr$^{-1}$): only the CRMHD galaxy maintains a steady outflow of $8 - 9$ M$_\odot$ yr$^{-1}$ while that in the MHD counterpart drops to $2$ M$_\odot$ yr$^{-1}$ and the HD one floats around $2 - 4$ M$_\odot$ yr$^{-1}$, with the odd jump to $8$ M$_\odot$ yr$^{-1}$ at $z \sim 2$. It is also important to note that in comparison with the HD and MHD runs, the
<figure><img src="image_10.png" /></figure>

**Figure 9.** Top panel: Outflow rates measured in a thin spherical shell located at a distance 0.2 $R_{vir,DM}$ from the galactic centre for the HD, MHD, and CRMHD runs. The histograms show the instantaneous mass-flow rate through the shell, while the solid lines represent the running median of these measurements using a filter of size 200 Myr. Middle panel: Star formation rate (SFR) of the galaxy averaged over time bins of 100 Myr. Bottom panel: Mass-loading factor $\eta$ computed using the running median of $\dot{M}$ and SFR shown in the top panels. The simulation with CRs has large outflow rates, × 10–100 higher than the HD and MHD simulations, and maintains them across cosmic time. Given the large suppression of star formation for $z \ge 3$, the CRMHD simulation also yields a mass-loading factor consistently exceeding unity, a value that is only reached in the HD and MHD simulations during the merger-induced starburst at $z \sim 2.3$.

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 a merger with a large companion, causing a significant increase in the cold gas mass of the galaxy (see Fig. 3) and a burst in star formation. We find that a corresponding large increase in the inflow rate on to the galaxy (Rodríguez Montero et al., in preparation), is fuelled by the cold, high density extraplanar gas supported by CR pressure that becomes unstable during the merger interaction. Overall, the cosmic evolution of $\dot{M}_{\text{outflow}}$ appears more steady than the SFR for the CRMHD run. This shows that while SN feedback is rarely efficient in the HD and MHD runs, even when large amounts of gas are converted into stars, in the CRMHD run a significant quantity of gas is consistently ejected from the galaxy even though fewer stars are formed overall.

<figure><img src="image_11.png" /></figure>

**Figure 10.** Histogram of mass outflow rates as a function of gas radial velocity measured in a thin shell located at 0.2 $R_{vir,DM}$ from the centre of the galaxy at $\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 snapshots over a $t_{dyn}$ period, using the method detailed in Section 3.2. Vertical dotted lines show the stack-weighted averages $v_{\text{esc}}$ for the redshift quoted, with the corresponding standard deviations denoted by the shaded regions. Values for the total mass outflow rate above $v_{\text{esc}}$ ($M_{\text{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. At high redshift, the shallower gravitational potential of the CRMHD galaxy allows a much larger fraction of the high velocity tail to exceed the local $v_{\text{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 M_{\odot} \text{yr}^{-1}$) occurring approximately at $z \sim 2.5$. The latter is characterized by a merger with a fairly massive companion. We compute stacked outflow rate histograms following the method detailed in Section 3.2 for weighted stacking of PDFs, employing snapshots within a dynamical time of the median redshift. These histograms show the total $\dot{M}$ separated into radial velocity bins for the same thin shell used in Fig. 9 located at a distance 0.2 $R_{vir,DM}$ from the galaxy centre. We show the HD, MHD, and CRMHD simulations in blue, purple, and green, respectively. The top panel of Fig. 10 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 the HD and MHD runs. We define the quantity $M_{\text{esc}}$ as the outflow rate of gas with a speed exceeding the local escape velocity. While the HD distribution is dominated by gas between 60 and 120 km s⁻¹, it rapidly drops at higher velocities, with just 0.01 M☉ yr⁻¹ exceeding the $v_{\text{esc}}$ of the galaxy (indicated by the vertical dotted blue line). Note that this instantaneous escape velocity is computed at each snapshot by counting the total mass enclosed within $0.2 R_{\text{vir,DM}}$. The MHD simulation instead shows a tail that does extend to higher velocities, contributing 10 times more to the mass outflow rate above $v_{\text{esc}}$ than HD. However, despite this higher efficiency in driving fast outflows, the majority of the outflowing mass between 60 km s⁻¹ and 120 km s⁻¹ remains the same, and the outflowing gas mass above $v_{\text{esc}}$ is still 3 times lower than in the CRMHD run. This is due to the fact that $v_{\text{esc}}$ is half the value for the CRMHD run than for the HD/MHD runs, because of the former's shallower potential well in the galactic region. This, in turn, is caused by the CR feedback persistently and significantly affecting the growth of the central stellar mass of NUT at higher redshifts. Therefore, the larger outflow rates in CRMHD at $z \sim 4$ can be explained by a larger fraction (almost 1 dex compared to HD and MHD) of 'low-velocity' (below 100 km s⁻¹) gas and an easier escape of fast outflows from a shallower gravitational potential. To further explore these trends, we also analyse the outflow velocity distribution at $z \sim 2.5$ (see bottom panel of Fig. 10). Although now the fast growth of NUT in the CRMHD simulation between $z = 3$ and $z = 2$ somewhat bridges the gap in $v_{\text{esc}}$, this simulation still has 6 times as much gas outflowing above $v_{\text{esc}}$ than the HD and MHD runs. This is not surprising given that during this period $\dot{M}_{\text{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_{\text{vir,DM}}$ lacks the kinetic energy necessary to escape the local gravitational potential in all simulations. However, it remains an open question as to whether this gas needs 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 systematically reaches higher mass loading factors than the no-CR runs. These outflows are mainly comprised of a slow moving wind which, without any further acceleration beyond $0.2 R_{\text{vir,DM}}$, will not be able to escape the halo and will eventually fall back onto the galaxy if left undisturbed by future feedback events.

### 3.3.3 Dynamics of multiphase outflows

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 disc. To provide evidence for this higher efficiency driving higher mass-loading factors in CRMHD, Fig. 11 shows the $\dot{M}_{\text{outflow}}$ and $\dot{M}_{\text{outflow}}$-weighted radial velocities of outflowing gas $v_{\text{outflow}}$ in the HD (blue), MHD (purple), and CRMHD (green) simulations, separated in different phases. HD and MHD outflows are dominated by the hot and ionized phases during the majority of the evolution of NUT. Outflowing gas in the warm neutral and cold phases in these runs are identified with abrupt, 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, in the CRMHD run, the hot phase never dominates the outflowing mass rate (see Fig. 12). It is the warm ionized phase which dominates $(1 - 10 M_{\odot} \text{yr}^{-1})$, with

considerable contributions of the warm neutral and cold gas at $z \ge 4$. We find that multiphase outflows ($\sim 1:2:4:2$ ratios for cold, warm neutral, warm ionized and hot phases, respectively) are only found consistently in the CRMHD simulation, and at high redshift ($z \ge 4$). The HD and MHD runs fail to ever produce such commensurable gas mass ratios across all phases.

To understand whether the outflowing gas is gravitationally bound, we compare its radial velocity with the local escape velocity at $0.2 R_{\text{vir,DM}}$ (indicated by dashed lines in the four relevant panels of Fig. 11). While a radial velocity higher than $v_{\text{esc}}$ is, strictly speaking, not a sufficient (or even necessary) condition for gas to completely escape the gravitational potential of a halo, it nevertheless constitutes a useful proxy, and is thus commonly used both in observational and theoretical studies. The radial velocities of the hot phase (top left-hand panel) range from $\sim 100 \text{ km s}^{-1}$ for $z \ge 3$ to $\sim 80 \text{ km s}^{-1}$ at lower redshifts in the simulations without CRs. For the case of CRMHD, the velocities are consistently higher, in the range $\sim 100 - 200 \text{ km s}^{-1}$. These outflow-weighted velocities are in agreement with observations of the warm ionized gas in local dwarf galaxies by McQuinn et al. (2019) and Marasco et al. (2023), while falling short of the $\sim 400 - 600 \text{ km s}^{-1}$ value reported in Concas et al. (2022) for their sample of galaxies between $1.2 < z < 2.6$.<fnref n="9" /> However, the CRMHD simulation is the only one consistently able to reach outflow velocities above the escape velocity within the galactic region, at least until $z \sim 2.5$, after which its velocity distribution becomes more similar to the non-CR cases. This behaviour is in agreement with the histograms shown in Fig. 10.

Warm neutral outflows are also observed in the HD and MHD simulations during merger periods, specifically during the active merger epoch between $6 \le z \le 8$ and around the late merger-induced starburst at $z \sim 2.3$. In Fig. 12, we show the fractional contribution of each phase to the total outflow rate through the previous thin shell located at $0.2 R_{\text{vir,DM}}$. While cold outflows are found in the HD and MHD runs, they only appear as rare, low velocity peaks in the $v_r$ panels. 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 outflows are linked to extended tidal features of satellites, rather than related to the rest of hot and warm ionized outflowing gas. In fact, CRMHD is the only run that consistently shows warm neutral and cold outflows across the majority of snapshots with non-negligible outflow velocities. These outflows are also spatially correlated with the hotter phases of the outflowing gas and are not dynamically linked to tidally stripped satellites (see Section 3.3.4 and Fig. 14 for further detail). Contrary to the hot and warm ionized phases, the outflow velocities of the warm neutral and cold gas evolve, increasing from $\sim 50 \text{ km s}^{-1}$ at high redshift to a maximum of $\sim 100 \text{ km s}^{-1}$ at $z \sim 2.5 - 2$. Similarly to the warm ionized phase, the median $v_r$ of these two phases sits well below $v_{\text{esc}}$ (except for peaks at $\sim 2.2$ and $\sim 2.5 \text{ Gyr}$ in the cold phase). Overall, we find that these cooler outflows struggle to reach escape velocities at $0.2 R_{\text{vir,DM}}$, even in the CRMHD simulation, despite its shallower gravitational potential.

<footnote n="9">⁹We note that outflows exceeding the escape velocity at a particular distance and time are neither a sufficient nor necessary condition for the eventual ejection of gas outside of the DM halo. Pressure forces in the outer halo can be larger than the inward gravitational pull, effectively accelerating the gas outwards, thus alleviating the need to reach the local escape velocity (e.g. Girichidis et al. 2018).</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Outflow budget and kinematics of hot (top left pair), warm ionized (top right pair), warm neutral (bottom left pair) and cold (bottom right pair) gas phases for the HD (blue), MHD (purple) and CRMHD (green) runs measured at 0.2 $R_{vir,DM}$. Each pair of panels shows the instantaneous outflow rate (top) and the $\dot{M}_{\text{outflow}}$-weighted radial velocity of the outflowing gas (bottom). The $\dot{M}_{\text{outflow}}$ panels show the same quantities as Fig. 9: the low opacity histogram is the instantaneous value, while the solid lines represent the running median over 200 Myr. In the $v_r$ panels, the escape velocity of each model at the shell radius is indicated by dashed lines, for ease of comparison with the median $v_{\text{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 ~100 km s⁻¹ for $z \ge 3$ and ~70 km s⁻¹ for $z \le 3$, dominating the outflow rate in both HD and MHD simulations. The hot phase in the CRMHD run reaches the largest outflow velocities (~150–200 km s⁻¹ for $z \ge 2.5$ and ~100 km s⁻¹ for $z \le 2.5$), with a non-negligible fraction of the outflowing gas mass exceeding $v_{\text{esc}}$ consistently up to $z \ge 2.5$. The dominant mass-flow rate comes from the warm ionized phase in the CRMHD simulation, while above $z \sim 5$, warm neutral and cold phases significantly contribute to the total outflow rate. By contrast, warm neutral and cold phases are scarce in the HD and MHD simulations, usually identified as low velocity peaks which match gas features from orbiting satellites crossing the shell. Note that, in the CRMHD run, while slower than the hot gas, these cold outflows can still reach ~100 km s⁻¹.</figcaption></figure>

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

The CRMHD simulation displays considerably higher outflow rates at 0.2 $R_{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_{vir,DM}$. This investigation

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_{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 in Fig. 13. For a pressured fluid $x$, this corresponds to $-\hat{r} \cdot \nabla P_i = \hat{r} \cdot (\rho \ddot{g})$, where $\hat{r}$ is the radial unit vector,
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Fractional contribution of the cold (blue lines), warm neutral (yellow lines), warm ionized (orange lines), and hot (red lines) phases to the total instantaneous outflow rates, for the HD (top panel), MHD (middle panel), and CRMHD (bottom panel) runs. The HD and MHD simulations are dominated by the hot and warm ionized phases for a similar amount of time. In contrast, in the CRMHD run, the warm ionized phase dominates the outflow rates almost at all times. The only exception is at very high redshifts (i.e. the first 0.5 Gyr of evolution), when the warm neutral and cold gas constitute a non-negligible fraction of the total outflow rates. All three simulations show an important contribution of warm neutral and cold outflows around $z \sim 6$, corresponding to an intense period of mergers for NUT.</figcaption></figure>

$\bar{g}$ is the local gravitational acceleration and $i$ 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_{\text{ther}}$, as solid circle symbols coloured in red linked by dashed lines, and CR pressure, $P_{\text{CR}}$, as green circles connected by solid lines. We describe in more detail how we compute pressure gradient support in Appendix B.

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 find that the hot and warm ionized phases show a large variance in these quantities, so we include the second and fourth quartiles of the underlying distribution as a shaded band in order to illustrate this. The three runs display a large amount of thermal pressure support against gravity for $z \ge 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 disc, 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 its CGM when compared to 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 simulations.<fnref n="10" /> 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. 11. While hot fast outflows in the HD and MHD simulations find it harder to travel through 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.

We compare in the right-hand panels of Fig. 13 the pressure support against gravity in the warm ionized, 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. This analysis confirms that the thermal pressure support (as well as the magnetic, which we do not include) in the HD and MHD runs is negligible and 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_{\text{vir,DM}}$. Due to outflow velocities lower than $v_{\text{esc}}$, outflows likely rain back down on to the galaxy as a galactic fountain. The CRMHD run features warm ionized outflows that are marginally supported against gravity, following an overall trend similar to that of the hot gas, but now with CR pressure gradients dominating over thermal pressure gradients at all times. In fact, CR pressure support dominates over thermal pressure support for the three phases in the right column of Fig. 13. The warm neutral and the cold components of the outflow show large accelerations due to CRs, exceeding by up to 7 times the gravitational acceleration. These phases display a large variability, also found in the analysis of $v_r$ (see Fig. 11). 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. 14, we show edge-on mass flow-weighted projections of NUT during a feedback event at $z = 2.41$ ($t \sim 2.9$ Gyr). This time was selected because of the presence of a corresponding peak in the outflow rate (i.e. $\sim 30 M_{\odot} \text{yr}^{-1}$, see Fig. 9) 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_{\text{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 $10^{-24} \text{ g cm}^{-3}$ in black (red) and at $10^{-22} \text{ g cm}^{-3}$ in red (blue). These panels uphold the canonical understanding of wind launching in disc galaxies: at the centre lies a cold and dense disc from which high velocity, hot outflows carve biconical perpendicular channels through the warm ionized CGM. This particular outflow is preceded by previous plumes of hot gas

<footnote n="10">¹⁰Nevertheless, 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.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Pressure gradient support against gravity for outflowing gas separated into different temperature phases for the 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 located at $0.2 R_{vir,DM}$ and divided into the phases described in Section 3.3.2. 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 ionized phase between $3 < z < 6$, but it supports and accelerates warm neutral and cold gas for $z \le 5$.</figcaption></figure>

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_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_{ther}$) and CR ($P_{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-hand panels), which display in both thermal and CR pressure support maps a clear biconical morphology with an opening angle of $\sim 30^\circ$. Dark red and blue colours represent high pressure acceleration antiparallel 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 disc 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. 13) as a method to identify galactic winds. In agreement with the peak observed at $t \sim 2.9$ Gyr in Fig. 13, 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$ 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 ionized phase of the gas is the dominant volume-filling phase of the outflow, and has a more spherical, 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
<figure><img src="image_15.png" /></figure>

**Figure 14.** Edge-on view of NUT in the CRMHD run during a strong feedback event at $z = 2.41$. Top panels: mass flow-weighted projections of a 20 kpc on a side cube centred on NUT showing the gas radial velocity (left-hand panel) and the gas temperature (right-hand panel). In the left-hand 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 disc, while warm gas at $\sim 10^{4.5}$ 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}$ g cm$^{-3}$ and $10^{-22}$ g cm$^{-3}$, respectively. The purple ring shows the 0.2 $R_{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 the text for detail). Each phase is plotted in two panels, corresponding to thermal (i.e. $-\nabla_r P_{ther}/(\rho g_r)$) and CR (i.e. $-\nabla_r P_{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 ionized (two top rightmost panels), warm neutral (two bottom leftmost panels), and cold (two bottom rightmost panels). To help visual analysis, we include in these panels the same gas density contours (green and blue) displayed on the temperature panel (top right-hand panel). The warm neutral and cold panels also include contour lines 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 outflowing waves. The warm ionized gas is the most volume filling phase, with a much more isotropic structure around NUT, and dominated by a uniform $-\nabla_r P_{CR}/(\rho g_r)$ profile. Warm neutral and cold outflows present a more clumpy and filamentary structure which quickly fades off at larger radii.
CR pressure. In fact, the warm ionized, warm neutral and cold phases exhibit the same support morphology, which suggests that they are all sourced by a common underlying CR pressurized halo that surrounds the galaxy. However, based on the mild CR pressure support of the warm ionized phase in Fig. 13, we claim that this phase will probably rain down as a galactic fountain given its low average radial velocity (see the top right-hand panel of Fig. 11). 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 disc, 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_{vir,DM}$ shell we employ to measure outflows (indicated by a dashed purple circle) in a time-scale shorter than our snapshot time resolution, naturally resulting in the fluctuating behaviour of warm neutral and cold phases present in Figs 11 and 13.

## 4 CONCLUSIONS

In this work, we have performed a detailed study of the effects of CR feedback on an MW-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 AMR, we simulated the multiphase ISM with a maximal spatial resolution of ~23 pc. Our HD simulation includes a model for a resolved ISM, with an MTT star formation prescription and a 'mechanical' SN feedback implementation of canonical strength. In our MHD simulation, we additionally followed the evolution of magnetic fields as they are amplified from a primordial seed of $3 \times 10^{-12}$ G. Our CRMHD simulation further accounted for CRs injected by SN explosions. 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 summarized as follows:

(i) Including CRs significantly reduces the resulting stellar mass of the galaxy 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 to boost SN energy deposition or to use calibrated galactic winds. This reduction in stellar mass is significant as it reaches approximately an order of magnitude at high redshift ($z \gtrsim 3$) before tapering to a factor ~4 at $z \sim 1.5$ (Section 3.1 and Fig. 4).

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

(iii) 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 environment where SN events

take place changes. In the CRMHD simulation, SN explosion sites exhibit a bi-modal distribution which both peaks at and reaches considerably lower densities than in the HD and MHD runs (see Fig. 7). This reduces the amount of SN energy lost due to cooling, increases the efficiency of this feedback mechanism and drives stronger outflows.

(iv) The mass outflow rates in the CRMHD simulation are significantly larger than those in the no-CR runs at all cosmic times, and yield a much higher mass loading factor $\eta$ especially at high redshift. In contrast to the simulations without CR where the dominant phase of outflowing gas mass is hot (i.e. $T > 10^5$ K), CRMHD outflows are dominated by the warm ionized phase (i.e. $9000 < T < 10^5$ K) and can even have a significant contribution from the cold phase (see Sections 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 multiphase galactic outflows expected soon.

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

We find CRs 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 individual roles is left to the next paper in this series, which will analyse the effects of CR feedback on the thermodynamic state of the CGM outside of the galactic region (Rodríguez Montero et al. in preparation).

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 (Katz et al. 2022; Nuñez-Castiñeyra et al. 2022) and different refinement criteria in the CGM (Bennett & Sijacki 2020; Rey et al. 2024). MW-like haloes 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, Sijacki & Smith 2022), or radiation from young stars (Rosdahl & Teyssier 2015; Emerick et al. 2018; Hopkins et al. 2020). 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 MW-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. SMA acknowledges support from the Kavli Institute for Particle Astrophysics and Cosmology
(KIPAC) Fellowship. 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 O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, ApJ, 770, 25

Agertz O. et al., 2021, MNRAS, 503, 5826

Armillotta L., Ostriker E. C., Jiang Y.-F., 2022, ApJ, 929, 170

Axford W. I., 1981, Ann. NY Acad. Sci., 375, 297

Baldry I. K., Sullivan T., Rani R., Turner S., 2020, MNRAS, 500, 1557

Barro G. et al., 2013, ApJ, 765, 104

Basu A., Roy S., 2013, MNRAS, 433, 1675

Beck R., 2016, A&AR, 24, 4

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 P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57

Bell E. F., 1978, MNRAS, 182, 147

Bennett J. S., Sijacki D., 2020, MNRAS, 499, 597

Biermann L., 1950, Z. Naturforsch, 5, 65

Binney J., Tremaine S., 2008, Galactic Dynamics, 2nd edn. Princeton Univ. Press, Princeton, p. 888

Blandford R., Eichler D., 1987, PhR, 154, 1

Blondin J. M., Wright E. B., Borkowski K. J., Reynolds S. P., 1998, ApJ, 500, 342

Booth C. M., Agertz O., Kravtsov A. V., Gnedin N. Y., 2013, ApJ, 777, L16

Boulares A., Cox D. P., 1990, ApJ, 365, 544

Bournaud F., Elmegreen B. G., Teyssier R., Block D. L., Puerari I., 2010, MNRAS, 409, 1088

Bournaud F. et al., 2013, ApJ, 780, 57

Breitschwerdt D., McKenzie J. F., Voelk H. J., 1991, A&A, 245, 79

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 M., Scannapieco E., 2020, ApJ, 905, 19

Buck T., Pfrommer C., Pakmor R., Grand R. J. J., Springel V., 2020, MNRAS, 497, 1712

Bustard C., Zweibel E. G., Bustard C., Zweibel E. G., 2021, ApJ, 913, 106

Camps P., Baes M., 2015, Astron. Comput., 9, 20

Cañas R., Elahi P. J., Welker C., Lagos C. d. P., Power C., Dubois Y., Pichon C., 2019, MNRAS, 482, 2039

Caprioli D., Spitkovsky A., 2014, ApJ, 794, 47

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 T. K., Kereš D., Gurvich A. B., Hopkins P. F., Trapp C., Ji S., Faucher-Giguère C.-A., 2022, MNRAS, 517, 597

Chisholm J., Tremonti Christy A., Leitherer C., Chen Y., 2016, MNRAS, 463, 541

Chisholm J., Tremonti C., Leitherer C., 2018, MNRAS, 481, 1690

Collins M. L., Read J. I., 2022, Nat. Astron., 6, 647

Commerçon B., Marcowith A., Dubois Y., 2019, A&A, 622, A143

Concas A. et al., 2022, MNRAS, 513, 2535

Crain R. A. et al., 2015, MNRAS, 450, 1937

Dashyan G., Dubois Y., 2020, A&A, 638, 123

Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827

de Blok W. J. G., 2010, Adv. Astron., 2010, 1

De Lucia G., Blaizot J., 2007, MNRAS, 375, 2

Dekel A., Birnboim Y., 2006, MNRAS, 368, 2

Dobbs C. L., Burkert A., Pringle J. E., 2011, MNRAS, 417, 1318

Dubois Y., Commerçon B., 2016, A&A, 585, A138

Dubois Y., Gavazzi R., Peirani S., Silk J., 2013, MNRAS, 433, 3297

Dubois Y., Peirani S., Pichon C., Devriendt J., Gavazzi R., Welker C., Volonteri M., 2016, MNRAS, 463, 3948

Dubois Y., Commerçon B., Marcowith A., Brahimi L., 2019, A&A, 631, A121

Dunkley J. et al., 2009, ApJS, 180, 306

Durrer R., Neronov A., 2013, A&AR, 21, 62

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 A., Bryan G. L., Low M.-M. M., 2018, ApJ, 865, L22

Ensslin T., 2004, J. Korean Astron. Soc., 37, 439

Evans II N. J., 1999, ARA&A, 37, 311

Evans N. J. et al., 2009, ApJS, 181, 321

Everett J. E., Zweibel E. G., Benjamin R. A., McCammon D., Rocks L., Gallagher III J. S., 2008, ApJ, 674, 258

Farber R., Ruszkowski M., Yang H.-Y. K., Zweibel E. G., 2018, ApJ, 856, 112

Farcy M., Rosdahl J., Dubois Y., Blaizot J., Martin-Alvarez S., 2022, MNRAS, 513, 5000

Farmer A. J., Goldreich P., 2003, ApJ, 604, 671

Federrath C., Klessen R. S., 2012, ApJ, 761, 156

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 G. B., Goldsmith D. W., Habing H. J., Field G. B., Goldsmith D. W., Habing H. J., 1969, ApJ, 155, L149

Fierlinger K. M. et al., 2016, MNRAS, 456, 710

Fromang S., Hennebelle P., Teyssier R., 2006, A&A, 457, 371

Geen S., Slyz A., Devriendt J., 2013, MNRAS, 429, 633

Geen S., Rosdahl J., Blaizot J., Devriendt J., Slyz A., 2015, MNRAS, 448, 3248

Gent F. A., 2012, Space Sci. Rev., 166, 281

Girichidis P. et al., 2016, ApJ, 816, L19

Girichidis P., Naab T., Hanasz M., Walch S., 2018, MNRAS, 479, 3042

Governato F. et al., 2010, Nature, 463, 203

Governato F. et al., 2012, MNRAS, 422, 1231

Guo F., Oh S. P., 2008, MNRAS, 384, 251

Haardt F., Madau P., 1996, ApJ, 461, 20

Hanasz M., Wóltański D., Kowalik K., 2009, ApJ, 706, L155

Hanasz M., Lesch H., Naab T., Gawryszczak A., Kowalik K., Wóltański D., 2013, ApJ, 777, L38

Hanasz M., Strong A. W., Girichidis P., 2021, Living Rev. Sol. Phys., 7, 2

Heckman T. M. et al., 1995, ApJ, 448, 98

Hennebelle P., Chabrier G., 2011, ApJ, 743, L29

Hennebelle P., Iffrig O., 2014, A&A, 570, A81

Hennebelle P. et al., 2022, A&A, 668, A147

Hollenstein L., Caprini C., Crittenden R., Maartens R., 2008, Phys. Rev. D, 77, 063517

Hopkins P. F., Quataert E., Murray N., 2011, MNRAS, 417, 950

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 P. F. et al., 2020, MNRAS, 492, 3465

Hopkins P. F., Chan T. K., Ji S., Hummels C. B., Kereš D., Quataert E., Faucher-Giguère C. A., 2021a, MNRAS, 501, 3640

Hopkins P. F., Squire J., Chan T. K., Quataert E., Ji S., Kereš D., Faucher-Giguère C.-A., 2021b, MNRAS, 501, 4184

Inoue S., Dekel A., Mandelker N., Ceverino D., Bournaud F., Primack J., 2016, MNRAS, 456, 2052

Ipavich F. M., 1975, ApJ, 196, 107

Ji S. et al., 2020, MNRAS, 496, 4221

Ji S. et al., 2021, MNRAS, 505, 259

Jóhannesson G., Porter T. A., Moskalenko I. V., 2019, ApJ, 879, 91

Jubelgas M., Springel V., Enßlin T., Pfrommer C., 2008, A&A, 481, 33
Katz H. et al., 2022, preprint (arXiv:2211.04626)

Kay S. T., Pearce F. R., Frenk C. S., Jenkins A., 2002, MNRAS, 330, 113

Kennicutt R. C., 1998, ARA&A, 36, 189

Kereš D., Katz N., Davé R., Fardal M., Weinberg D. H., 2009, MNRAS, 396, 2332

Kim C. G., Ostriker E. C., 2015, ApJ, 802, 99

Kim C.-G., Ostriker E. C., 2018, ApJ, 853, 173

Kimm T., Cen R., 2014, ApJ, 788, 121

Kimm T., Cen R., Devriendt J., Dubois Y., Slyz A., 2015, MNRAS, 451, 2900

Kimm T., Katz H., Haehnelt M., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS, 466, 4826

Kimm T., Haehnelt M., Blaizot J., Katz H., Michel-Dansac L., Garel T., Rosdahl J., Teyssier R., 2018, MNRAS, 475, 4617

Koudmani S., Sijacki D., Smith M. C., 2022, MNRAS, 516, 2112

Kroupa P., 2001, MNRAS, 322, 231

Krumholz M. R., Tan J. C., 2007, ApJ, 654, 304

Kulpa-Dybel K. et al., 2015, A&A, 575, A93

Kulsrud R., Pearce W. P., 1969, ApJ, 156, 445

Kuwabara T., Ko C.-M., 2020, ApJ, 899, 72

Leitherer C. et al., 1999, ApJS, 123, 3

Leroy A. K. et al., 2023, ApJ, 944, L10

Lopez-Rodriguez E., Guerra J. A., Asgari-Targhi M., Schmelz J. T., 2021, ApJ, 914, 24

Maiolino R., Mannucci F., 2019, A&AR, 27, 3

Mannucci F., Cresci G., Maiolino R., Marconi A., Gnerucci A., 2010, MNRAS, 408, 2115

Mao S. A., Ostriker E. C., 2018, ApJ, 854, 89

Marasco A. et al., 2023, A&A, 670, A92

Martin-Alvarez S., Devriendt J., Slyz A., Teyssier R., 2018, MNRAS, 479, 3343

Martin-Alvarez S., Slyz A., Devriendt J., Gómez-Guijarro C., 2020, MNRAS, 495, 4475

Martin-Alvarez S., Katz H., Sijacki D., Devriendt J., Slyz A., 2021, MNRAS, 504, 2517

Martin-Alvarez S., Sijacki D., Haehnelt M. G., Farcy M., Dubois Y., Belokurov V., Rosdahl J., Lopez-Rodriguez E., 2023, MNRAS, 525, 3806

Mashchenko S., Couchman H. M. P., Wadsley J., 2006, Nature, 442, 539

McQuinn K. B. W., van Zee L., Skillman E. D., 2019, ApJ, 886, 74

Miyoshi T., Kusano K., 2005, J. Comput. Phys., 208, 315

Mocz P., Pakmor R., Springel V., Vogelsberger M., Marinacci F., Hernquist L., 2016, MNRAS, 463, 477

Morlino G., Caprioli D., 2012, A&A, 538, A81

Moster B. P., Naab T., White S. D. M., Moster B. P., Naab T., White S. D. M., 2013, MNRAS, 428, 3121

Mulcahy D. D. et al., 2014, A&A, 568, A74

Murray N., Quataert E., Thompson T. A., 2005, ApJ, 618, 569

Naab T., Ostriker J. P., 2017, ARA&A, 55, 59

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 A., Grenier I. A., Bournaud F., Dubois Y., Youssef F. R. K., Hennebelle P., 2022, preprint (arXiv:2205.08163)

Pakmor R., Springel V., 2013, MNRAS, 432, 176

Pakmor R., Pfrommer C., Simpson C. M., Springel V., 2016, ApJ, 824, L30

Peschken N., Hanasz M., Naab T., Wóltański D., Gawryszczak A., 2023, MNRAS, 522, 5529

Pillepich A. et al., 2018, MNRAS, 473, 4077

Powell L. C., Slyz A., Devriendt J., 2011, MNRAS, 414, 3671

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 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 R. E., Silk J., Pudritz R. E., Silk J., 1989, ApJ, 342, 650

Quataert E., Jiang Y.-F., Thompson T. A., 2021a, MNRAS, 510, 920

Quataert E., Thompson T. A., Jiang Y.-F., 2021b, MNRAS, 510, 1184

Rasera Y., Teyssier R., 2006, A&A, 445, 1

Rey M. P., Katz H. B., Cameron A. J., Devriendt J., Slyz A., 2024, MNRAS, 528, 5412

Rhee J., Elahi P., Yi S. K., 2022, ApJ, 927, 129

Rieder M., Teyssier R., 2017, MNRAS, 471, 2674

Rodríguez Montero F., Martin-Alvarez S., Sijacki D., Slyz A., Devriendt J., Dubois Y., 2022, MNRAS, 511, 1247

Rosdahl J., Teyssier R., 2015, MNRAS, 449, 4380

Rosdahl J. et al., 2018, MNRAS, 479, 994

Rosen A., Bregman J. N., 1995, ApJ, 440, 634

Ruszkowski M., Yang H.-Y. K., Zweibel E., 2017, ApJ, 834, 208

Salem M., Bryan G. L., Hummels C., 2014, ApJ, 797, L18

Schlickeiser R., Schlickeiser R., 2009, MNRAS, 398, 1483

Schmidt M., 1959, ApJ, 129, 243

Semenov V. A., Kravtsov A. V., Caprioli D., 2021, ApJ, 910, 126

Shadmehri M., 2009, MNRAS, 397, 1521

Siejkowski H., Soida M., Otmianowska-Mazur K., Hanasz M., Bomans D. J., 2009, A&A, 510, A97

Silk J., 2010, Proc. Int. Astron. Union, 6, 273

Skilling J., 1971, ApJ, 170, 265

Skilling J., 1975, MNRAS, 172, 557

Smith M. C., Sijacki D., Shen S., 2018, MNRAS, 478, 302

Smith M. C., Sijacki D., Shen S., 2019, MNRAS, 485, 3317

Socrates A., Davis S. W., Ramirez-Ruiz E., 2006, ApJ, 687, 202

Somerville R. S., Davé R., 2015, ARA&A, 53, 51

Springel V., Hernquist L., 2003, MNRAS, 339, 312

Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn T., 2006, MNRAS, 373, 1074

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 L. K., Stevens I. R., Strickland D. K., Heckman T. M., 2003, MNRAS, 342, 690

Tasker E. J., 2011, ApJ, 730, 11

Teyssier R., 2002, A&A, 385, 337

Teyssier R., Fromang S., Dormy E., 2006, J. Comput. Phys., 218, 44

Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MNRAS, 429, 3068

Thilker D. A. et al., 2023, ApJ, 944, L13

Thomas T., Pfrommer C., Pakmor R., 2023, MNRAS, 521, 3023

Thornton K., Gaudlitz M., Janka H., Steinmetz M., 1998, ApJ, 500, 95

Tillson H., Devriendt J., Slyz A., Miller L., Pichon C., 2015, MNRAS, 449, 4363

Tortora C., Hunt L. K., Ginolfi M., 2022, A&A, 657, A19

Tóth G., Tóth Gábor, 2000, JCoPh, 161, 605

Tremonti C. A. et al., 2004, ApJ, 613, 898

Tweed D., Devriendt J., Blaizot J., Colombi S., Slyz A., 2009, A&A, 506, 647

Vachaspati T., 1991, Phys. Lett. B, 265, 258

Vachaspati T., 2021, Rep. Prog. Phys., 84, 074901

Vazfa F., Jones T. W., Brüggen M., Brunetti G., Gheller C., Porter D., Ryu D., 2017, MNRAS, 464, 210

Veilleux S., Cecil G., Bland-Hawthorn J., 2005, ARA&A, 43, 769

Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nat. Rev. Phys., 2, 42

Wadepuhl M., Springel V., 2011, MNRAS, 410, 1975

Wagner A. Y., Falle S. A. E. G., Hartquist T. W., Pittard J. M., 2005, A&A, 430, 567

Werhahn M., Pfrommer C., Girichidis P., 2021, MNRAS, 508, 4072

Whalen D., van Veelen B., O'Shea B. W., Norman M. L., 2008, ApJ, 682, 49

Williams J. P., McKee C. F., Williams J. P., McKee C. F., 1997, ApJ, 476, 166

Wise J. H., Turk M. J., Norman M. L., Abel T., 2012, ApJ, 745, 50

Wolfire M. G. et al., 1995, ApJ, 443, 152

Xu S., Lazarian A., 2022, ApJ, 927, 94

Zubko V., Dwek E., Arendt R. G., 2004, ApJS, 152, 211

Zuckerman B., Evans N. J. I., 1974, ApJ, 192, L149

Zweibel E. G., 2017, Phys. Plasmas, 24, 055402

## APPENDIX A: RADIAL PROFILES OF PRESSURE COMPONENTS

In this appendix, we further review the dominance of the energy components explored in this work through volume-weighted profiles.
<figure><img src="image_16.png" /><figcaption>**Figure A1.** Volume-weighted radial profiles of the thermal (red), magnetic (purple) and CR (green) pressures centred on NUT at $z = 1.5$ for the CRMHD simulation. CR pressure is the dominant component from the inner galactic disc ($r \sim 0.15$ kpc) to the galaxy-halo interface ($r \sim 20$ kpc). Thermal pressure only dominates within the very high density gas at the centre of the galaxy, where magnetic and CR pressure approach equipartition. Magnetic pressure decreases more rapidly than thermal or CR pressures, thus remaining subdominant at all outer radii. The vertical purple region indicates the location of the thin shell at $0.2 R_{\text{vir,DM}}$. We also include the volume-weighted gas density profile (dashed black line) for context.</figcaption></figure>

This is done at $z = 1.5$ (see projections in Fig. 2), when NUT is well settled as a disc galaxy and no recent galaxy mergers have taken place. In Fig. A1, we show the radial profiles of thermal (red), magnetic (purple), and CR (green) pressures extending from the centre of the galaxy to slightly beyond $0.2 R_{\text{vir,DM}}$. These quantities are volume-weighted to convey the properties of the volume-filling gas above the disc and at the galaxy-halo interface. To aid in the interpretation across these scales, we have also included the volume-weighted radial density profile. As expected for such weighting, this profile follows a single power-law trend up to $r \sim 10$ kpc, where a knee can be attributed to the edge of the cold disc (Fig. 2). CR pressure dominates over the other pressure components within the majority of the galactic region. Only inside the innermost central region of NUT does thermal pressure exceed CR pressure by about a factor of 4. This is evidence for a CR dominated halo, capable of supporting gas in hydrostatic equilibrium. In order to connect these results with the thin shell analysis performed throughout this work, we have included a purple vertical bar which represents the fiducial shell used for outflow measurements, which for this particular redshift is located at $\sim 20$ kpc. Within this shell, all pressure profiles display negative gradients, suggesting outward acceleration against the local gravitational pull. Properly studying this claim requires careful examination of the pressure balances on a cell-per-cell basis, as we have done with our thin shell analysis (see Appendix B). Magnetic pressure is subdominant across virtually all radii, falling off with a radial trend similar to that of the gas density. At distances larger than the galactic disc radius ($r \gtrsim 10$ kpc), magnetic pressure decreases even more rapidly, rendering it of secondary importance for outflow dynamics close to $0.2 R_{\text{vir,DM}}$.

## APPENDIX B: THIN SHELL ANALYSIS

In the present work, we employ spherical thin shells to investigate the properties of outflows. This simple method is a powerful tool to study the role of CR and particularly effective in an explicit cosmological context such as that of our NUT simulations. This appendix is focused

on further describing the subtleties of our method, and on putting it in perspective with the more traditional approach of radial profiles used in the literature (e.g. Hanasz et al. 2013; Girichidis et al. 2016, 2018; Semenov, Kravtsov & Caprioli 2021; Thomas et al. 2023).

### B1 Shell width convergence study

We take a purely Eulerian approach for our measurements of outflow properties. To preserve the distinction between the galactic region and the outer halo of NUT, we choose the galaxy-halo interface as our physically motivated fixed spatial location to perform our outflow measurements. Outflow rates are computed using equation (9), which requires a choice for the shell width $\delta r_{\text{shell}}$. Choosing the most appropriate value of $\delta r_{\text{shell}}$ entails balancing two competing considerations: (1) equation (9) is derived from infinitesimal mass conservation, hence the measurement becomes more accurate for thinner shells, however (2) with a thicker shell the shot noise observed for discontinuous gas flows (e.g. warm neutral and cold phases in Fig. 14) can be reduced. We review this balance in Fig. B1 by studying the influence of $\delta r_{\text{shell}}$ on $\dot{M}_{\text{outflow}}$ and $-\nabla_r P_i / \rho g_r$ for the four phases explored in Fig. 13. We compare outflow related quantities for measurements performed with the fiducial shell width $\delta r_{\text{shell}} = 0.02 R_{\text{vir,DM}}$ (green lines) to those performed with a shell two times wider (dark olive lines) and a shell four times thinner (dark golden lines). The top panel of Fig. B1 shows that the outflow rates are well converged over a large range of shell widths around our fiducial $\delta r_{\text{shell}}$. In the next 4 panels, we reproduce the results for the support of the CR pressure gradient against gravity studied in Fig. 13. As predicted, the most volume-filling phases (hot and warm ionized) show no difference in the averaged pressure support against gravity for different shell widths. However, for the warm neutral and cold phases, the shell thickness has a non-negligible influence on the maxima and minima measured. Since the thinnest shell (i.e. $\delta r_{\text{shell}} = 0.005 R_{\text{vir,DM}}$) encompasses a significantly smaller volume, the chances of capturing cold ‘clumps’ crossing the shell are drastically reduced. Hence, a thinner shell leads to a greater stochasticity in the measured support, although the average temporal trend is still consistent across measurements with different shell width.

### B2 Comparing the shell method to radial profiles

To compute the ratio $-\nabla_r P_i / \rho g_r$, we need to compute pressure gradients and gravitational acceleration. The gravitational acceleration at each cell centre is directly provided as an output from RAMSES. Pressure gradients are computed in post-processing as they would be computed on-the-fly by RAMSES. More specifically, local pressure gradients are computed by walking down the AMR tree until leaf cells (i.e. cells which are not further refined) are reached. The gradient of quantity X for cell j in the x-direction is given by

$$ \partial_x X_j = \frac{X_{j+1/2} - X_{j-1/2}}{\Delta x_{j+1/2} + \Delta x_{j-1/2}}, \quad (B1) $$

where $j + 1/2$ and $j - 1/2$ indicate the left and right neighbouring cells along x, respectively. The quantities in the denominator $\Delta x_{j+1/2}$ and $\Delta x_{j-1/2}$ provide the distance in the x direction from the centres of the neighbouring cells to the centre of cell j. If all neighbouring cells are on the same level l as the cell j, the distance is simply the cell size for level l, $\Delta x^l$. However, if one of the neighbours is a cell at a lower level l - 1, its cell size is twice as large, hence the distance along x would instead be $1.5 \Delta x^l$. Following this procedure also along the y and z directions, the gradient $\nabla X$ is obtained for each cell. By
<figure><img src="image_17.png" /></figure>

**Figure B1.** Exploration of the influence of shell thickness $\delta r_{\text{shell}}$ on the outflow-related quantities studied in this work. The adopted fiducial width $\delta r_{\text{shell}} = 0.02 R_{\text{vir,DM}}$ is represented by the green lines, while twice as thick (dark olive lines) and four times thinner (dark golden lines) are also shown for comparison. All shells are situated at $0.2 R_{\text{vir,DM}}$. Top panel: Outflow rate, with histograms showing the instantaneous mass-flow rate through the thin shell, while the solid lines are the running median calculated over a 200 Myr period. Bottom panels: CR pressure gradient support against gravity for outflowing gas, separated into the outflow phases explored in this work (from top to bottom): hot, warm ionized, warm neutral, and cold. In the shaded region, gravity dominates the dynamics, while in the unshaded region, gas is accelerated outwards by $\nabla P_{\text{CR}}$. Outflow rates are well converged for different shell widths, as is the pressure gradient support for the volume-filling phases (hot and warm ionized). A thinner shell gives rise to higher stochasticity in the measured support of warm neutral and cold phases.

computing local gradients in this manner, we maintain consistency with the numerical CRMHD implementation.

To explore whether our thin shell method yields results consistent with those obtained from profiles, in Fig. B2 we present an edge-on view of the outflowing gas during the burst event at $z \sim 2.41$ (the same redshift as in Fig. 14). In this case, the galactic angular momentum $j_{\text{gal}}$ is directed to the right of the page, thus the region to the left of the vertical dashed purple line is below the galactic plane, while the region to the right is above it. All quantities (both projected and binned along the profile) are average $\dot{M}_{\text{outflow}}$-weighted,

<figure><img src="image_18.png" /></figure>

**Figure B2.** Thin shell method comparison vis-à-vis projections and radial profiles. Top panels: $\dot{M}$-weighted, 6 kpc thin projections of the outflowing gas in NUT during the burst event at $z \sim 2.41$. In all projections, the galaxy is observed edge-on, with the galactic angular momentum $j_{\text{gal}}$ pointing towards the right edge of the figure. In the top two panels, this side-on view reveals fast, hot, bipolar outflows emanating from the galactic disc (dashed purple vertical line). In the subsequent two panels we display only half of the galactic region for each phase (black dashed line is the central axis). These projections show the support against gravity of the $P_{\text{CR}}$ gradient for hot, warm ionized, warm neutral and cold phases. We have overlaid white contours where $-\nabla_r P_{\text{CR}} / \rho g_r \ge 1$. Green and red contours indicate the projected gas density in the galactic disc as in Fig. 14. The purple ring shows the $0.2 R_{\text{vir,DM}}$ thin shell. Bottom panels: Radial profiles of $P_{\text{ther}}$ (dashed) and $P_{\text{CR}}$ (solid) gradients support for the top (right-hand side bottom panel) and the bottom (left-hand side bottom panel) galactic hemispheres. We distinguish the hot (red), warm ionized (orange), warm neutral (yellow), and cold (blue) gas phases as in the top projections. Symbols show the resulting support in each gas phase as measured by the thin shell method, with the filling colour showing the type of pressure gradient support (red for thermal, green for CRs) and the colours of the outlines of the symbols and error bars indicating the gas phase. Error bars show the upper and lower quartiles for the underlying distribution within the thin shell. The shell method is able to capture the $P_{\text{CR}}$ gradient support of hot, warm neutral and cold outflows, and the insufficient support for the warm ionized gas.
i.e. they only are presented for gas that is outflowing as seen in the reference frame of the galaxy. The two top panels show the radial velocity and the gas temperature, with fast outflows clearly seen as bi-conical hot regions emanating from the galactic disc. The slowly (~40–50 km/s) outflowing gas is mainly co-spatial with the colder gas at ~10<sup>4.5</sup> K. In both panels we have overlaid gas density contours as in Fig. 14 (in red and green), to delineate the extent of the galactic cold disc. Progressing downwards in the figure, we show the support against gravity provided by −∇<sub>r</sub>P<sub>CR</sub> in the hot, warm ionized, warm neutral and cold phases. Given the significant axisymmetry observed in the radial velocity and temperature projections (see the symmetry about the horizontal white dashed line in upper two panels), in the third and fourth panels we have opted to show the phase-separated projections for only half of the galactic region, with the symmetry axis, represented by the horizontal black dashed line, dividing the galactic region in two. Different from Fig. 14, here we only show cells for which the acceleration by the P<sub>CR</sub> gradient works against the radially inward gravitational acceleration. Additionally, we have added white contours in these projections to show where ∇<sub>r</sub>P<sub>CR</sub> is equal or larger than ρ<sub>g</sub>r. White regions in the projections indicate cells for which none of the conditions stated above are met.

In the bottom panel of Fig. B2, we present the summary statistics that quantify the features and trends observed in the projections. We have obtained spherical radial profiles for gas within the two hemispheres above and below the galactic plane, out to 0.2*R<sub>vir,DM</sub>. For each of the outflowing gas phases explored in this work (see Section 3.3.3), namely hot (red lines), warm ionized (orange lines), warm neutral (yellow lines), and cold (blue lines), we present radial profiles of thermal (dashed lines) and CR (solid lines) pressure gradient support against gravity. Hot gas profiles show that both thermal and CR pressure gradients are capable of accelerating hot gas, but thermal pressure gradients are more important close to the galactic disc while CR pressure gradients provide a more consistent acceleration across all distances from the galaxy. This hints at the underlying acceleration mechanisms driven by the two pressure components: thermal acceleration is due to large pressure contrasts between the volume filling warm gas and the injected hot outflowing gas, while CR acceleration is caused by a quasi-uniform CR halo permeating the galactic region. Corroborating the lack of white contours in the warm ionized projection, neither thermal nor CR pressure gradient radial profiles for the warm ionized phase show efficient radial acceleration against gravity (i.e. −∇<sub>r</sub>P<sub>i</sub>/ρ<sub>g</sub>r ≤ 1). Warm neutral and

cold gas profiles show that cool ‘cloudlets’ become more efficiently
accelerated away from the inner region of the galactic region.

We compare the results of these radial profiles to the averaged
measurements in our thin shell at 0.2 R<sub>vir,DM</sub>. The spatial scale
of this shell, with our fiducial thickness of 0.02R<sub>vir,DM</sub>, is shown
in the projections of Fig. B2 as a low-opacity purple ring. Shell-
averaged values are indicated, as filled circle symbols with error
bars, in the lower panel together with the radial profiles. We caution
that the symbols record values obtained by averaging over the
whole shell at ~ 15.6 kpc. To avoid overcrowding the panel, the
CR pressure symbols (green-filled markers) are only shown on the
plot of the radial profiles for the hemisphere above the galactic plane
(right bottom panel), while the thermal pressure symbols (red-filled
markers) are only shown in the plot of the radial profiles for the lower
hemisphere (left bottom panel). In agreement with the projections and
the profiles, the shell method attributes negligible thermal support
for all outflowing phases, except for the hot phase. For the case
of CR pressure gradients, the shell method again agrees with the
profiles in indicating CR support for the hot, warm neutral and
cold phases, while -∇<sub>r</sub>P<sub>CR</sub>/ρg<sub>r</sub> for the warm ionized phase remains
with an average value below unity in both hemispheres. The radial
profiles for the warm neutral and cold phases display strong jumps in
their pressure gradient support. Unlike for the hot or warm ionized
volume-filling phases, jumps in the warm neutral and cold profiles
are caused by ‘clumps’ in these outflows. In fact, it is easy to match
these features in the profiles to discontinuities in the projections
for these phases. Therefore, if these discontinuities move from their
launching point near the disc to the outer radius of the galactic region,
time-averaged radial profiles will not be able to detect these time-
varying features. While the thin shell method cannot fully grasp the
complexity of individual outflow profiles, it is capable of capturing
the outflow physical behaviour. The inhomogeneous nature of cooler
gas phases and the varying star formation history across the extended
cosmological evolution of the NUT galaxy justify and explains the
strength of our thin shell method over profiles traditionally used
to study established, smoothly star-forming, isolated galaxies. By
maintaining a spherical geometry for the shell across all cosmic time,
we maximize the adaptability of our method to galaxy geometry,
outflow direction and gas phase.

This paper has been typeset from a TeX/LATeX file prepared by the author.