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

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


INTRODUCTION
Realistically modelling galaxy formation within the ΛCDM cosmological paradigm has proven to be an extremely complex task, as its multi-scale and multi-physics nature quickly develops a non-linear behaviour that is difficult to grasp and model.The most sophisticated amongst modern cosmological simulations have confronted this challenge by explicitly solving the equations governing the evolution of matter across the vast range of scales spanning from the cosmic large-scale structure down to the interstellar medium (ISM) of individual galaxies.For processes such as star formation which are unresolved, yet indispensable for galaxy formation, numerical prescriptions are implemented.Much of their success (for a review see e.g.Vogelsberger et al. 2020, and references therein) resides in their recognition that, without the injection of energy by clustered supernova (SN) events, the gravitational collapse of cold-dense gas leads to the consumption of the gas reservoir over a free-fall time, resulting in too-efficient star formation (Kay et al. 2002; Bournaud ★ E-mail: currodri@gmail.comet al. 2010;Dobbs et al. 2011;Hopkins et al. 2011;Tasker 2011;Silk 2010;Moster et al. 2013) in stark difference with observations (e.g.Zuckerman & Evans 1974;Williams et al. 1997;Kennicutt 1998;Evans II 1999;Krumholz & Tan 2007;Evans et al. 2009).Besides the need for reproducing the inefficient formation of stars (e.g.Stinson et al. 2006;Teyssier et al. 2013;Springel & Hernquist 2003), including such a considerable energy injection through SN or active galactic nuclei (AGN) establishes a feedback loop which is of utter importance for reproducing many global observables of the galaxy population.Illustrative examples are the mass-metallicity relation (Tremonti et al. 2004;Kereš et al. 2009;Mannucci et al. 2010;Maiolino & Mannucci 2019;Tortora et al. 2022), the metal enrichment of the circum-(CGM) and intergalactic medium (IGM) (e.g.Veilleux et al. 2005;Brooks et al. 2007;Chisholm et al. 2018), the quenching and morphological transformation of late-type galaxies into early-type galaxies (e.g.De Lucia & Blaizot 2007;Naab et al. 2007;Dubois et al. 2013Dubois et al. , 2016)), the flattening of inner dark matter profiles (e.g.Mashchenko et al. 2006;Governato et al. 2010;de Blok 2010;Teyssier et al. 2013) and the formation of bulge-less galaxies by expelling low angular momentum gas (Governato et al. 2012).
However, this success is built on a phenomenological strategy, in which subgrid models for feedback processes are calibrated against local observed scaling relations (e.g.Crain et al. 2015;Pillepich et al. 2018;Rosdahl et al. 2018;Davé et al. 2019) thus considerably curtailing predictive power.For stellar feedback, state-of-the-art simulations are now moving into a regime where the multi-phase structure of the interstellar medium is beginning to be spatially resolved.This has driven the development of more physically motivated treatments of SNe feedback (e.g.Hopkins et al. 2014;Kimm & Cen 2014;Kimm et al. 2015), as well as other forms of stellar feedback that could not be accounted for in the previous generations of numerical simulations such as e.g.stellar radiation, including photoionisation heating and radiation pressure (Murray et al. 2005;Hopkins et al. 2011;Agertz et al. 2013;Rosdahl & Teyssier 2015;Emerick et al. 2018;Smith et al. 2018Smith et al. , 2019)), stellar winds (Hopkins et al. 2011;Agertz et al. 2013;Geen et al. 2015;Fierlinger et al. 2016;Smith et al. 2019), or direct pressure from Lyman- photons (Kimm et al. 2018).Although simulations are beginning to resolve the scales and physical processes required to model star formation through prescriptions aimed to reproduce their formation in giant molecular clouds, almost all aspects of the actual launching of outflows at galactic scales and their thermodynamic properties remains a source of intense debate in the theoretical and numerical simulations community.Furthermore, the observations that are used to constrain stellar feedback models (e.g.Collins & Read 2022;Leroy et al. 2023;Thilker et al. 2023) suffer from large limitations, with reported mass-loading factors (defined as the ratio of gas mass ejected from the galaxy to the current rate of gas depletion by star formation) varying by more than 2 orders of magnitude in galaxies of similar mass (e.g.Veilleux et al. 2005;Chisholm et al. 2016).There are also considerable difficulties involved in detecting the emission from hot (∼ 10 6 K) gas (Heckman et al. 1995;Summers et al. 2003;McQuinn et al. 2019;Marasco et al. 2023;Concas et al. 2022), which theoretical models argue carry a significant fraction of the kinetic energy of the wind (e.g.Kim & Ostriker 2018).
Cosmic rays (CRs), relativistic particles which are accelerated by diffusive shock acceleration (e.g.Axford 1981;Bell 1978;Blandford & Eichler 1987), act as a non-thermal source of energy, and pioneering numerical investigations (Jubelgas et al. 2008;Wadepuhl & Springel 2011;Booth et al. 2013;Salem et al. 2014, e.g.) of their impact indicate that they can have a significant impact on star formation.Observations of the Milky Way and nearby galaxies suggest that CRs and magnetic fields are in equipartition with the thermal and turbulent energies in their mid-plane (Boulares & Cox 1990;Basu & Roy 2013;Beck 2016;Zweibel 2017;Lopez-Rodriguez et al. 2021).Furthermore, the peak of their spectrum at ∼ 1 GeV leads them to exhibit smooth spatial gradients on ≳ 1 kpc scale lengths (e.g.Ensslin 2004).Consequently, for some time now they have been put forward as a fundamental player in the dynamics of extra-planar gas and the launching of outflows (e.g.Ipavich 1975;Breitschwerdt et al. 1991;Ptuskin et al. 1997;Socrates et al. 2006;Everett et al. 2008;Mao & Ostriker 2018;Quataert et al. 2021b,a).These ∼ 1 GeV CRs behave as a relativistic fluid with a soft equation of state (i.e.adiabatic index  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 et al. 2009) than the radiative cooling of thermal gas.They rapidly diffuse away from their injection sites due to diffusion lengths larger than typical molecular clouds (∼ 100 pc), reaching the more dilute gas immediately above the galactic disk, which is easier to accel-erate against the gravitational potential.Therefore, in the past few years there has been an increased effort in modelling CR feedback in numerical studies of galactic outflows, showing broadly consistent trends across hydrodynamic schemes and numerical implementations of the CR Fokker-Planck equation (Skilling 1971;Hanasz et al. 2021, and references therein): (i) CRs thicken the gaseous disk by altering the pressure-gravity balance, changing the morphology of the galaxy (Booth et al. 2013;Salem et al. 2014;Ruszkowski et al. 2017;Dashyan & Dubois 2020;Chan et al. 2022;Farcy et al. 2022;Nuñez-Castiñeyra et al. 2022;Buck et al. 2020;Martin-Alvarez et al. 2022).In addition, CR feedback significantly affects the star formation history of galaxies by (ii) reducing star formation rates (e.g.Jubelgas et al. 2008;Booth et al. 2013;Pakmor et al. 2016;Chan et al. 2019;Dashyan & Dubois 2020;Hopkins et al. 2020) and by (iii) launching more heavily mass-loaded outflows, these being denser and colder (e.g.Booth et al. 2013;Hanasz et al. 2013;Hopkins et al. 2021a;Farcy et al. 2022;Peschken et al. 2022).
This modelling progress has established that CRs in the ∼ 1 GeV regime are likely an important ingredient of galaxy formation.However, a clear understanding of its coupling with other physical processes of relevance remains far from being explored.While isolated galaxy simulations find CRs modelled with a constant diffusion coefficient  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 CRdriven outflows and cosmic gas inflows in a realistic CGM, that can be captured only by fully self-consistent cosmological simulations (see e.g.Ji et al. 2020Ji et al. , 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 MHD turbulence (e.g.Farmer & Goldreich 2003;Xu & Lazarian 2022), giving rise to a variation of up to 5 orders of magnitude in the CR mean free path (Armillotta et al. 2022) and to changes in the acceleration of cold clouds in the CGM (Brüggen & Scannapieco 2020;Bustard et al. 2021).Moreover, CRs can affect the evolution of magneto-thermal instabilities in galactic discs depending on the transport prescription and diffusion constant (e.g.Wagner et al. 2005;Shadmehri 2009;Kuwabara & Ko 2020), strongly altering its structure (Commerçon et al. 2019;Dashyan & Dubois 2020;Nuñez-Castiñeyra et al. 2022;Martin-Alvarez et al. 2022) and the conditions in which star formation and stellar feedback take place.Hence, an accurate modelling of the ISM is fundamental for a better understanding of CR feedback together with the launching and acceleration of CR-driven multiphase outflows (Girichidis et al. 2016;Kim & Ostriker 2018;Salem et al. 2014;Hopkins et al. 2021a;Martin-Alvarez et al. 2022), which remains to be carefully studied in a cosmological context.This is the first paper in a series in which we study the effects of CR feedback in galaxies with a resolved multi-phase ISM using highresolution cosmological zoom-in simulations, with the aim to further our understanding of the role CRs play in shaping galaxy properties.In Section 2 we outline our Milky Way-like galaxy (called nut) formation simulation set-up, and our CR and magnetohydrodynamics modelling.We then discuss the global effects of magnetic fields and CRs in the evolution of nut in Section 3.1, demonstrating how the CRs reduce the stellar mass by decreasing the star formation efficiency in Section 3.2 and increasing the SN feedback efficiency in Section 3.3 via modifications to the outflow phases and their more efficient acceleration.In Section 4 we summarise our main results, discuss their limitations and the prospect for future studies of CR feedback.

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

Cosmic ray magneto-hydrodynamics
We make use of the implementation of CRs from Dubois & Commerçon (2016) and Dubois et al. (2019), in which CRs are treated as a separate fluid represented by one additional energy equation in the set of MHD fluid equations: where  is the gas mass density, ì  is the gas velocity vector, ì  the gravitational field, ì is the streaming velocity (where ì  A is the Alfvén velocity and ì  = ì /| ì | the magnetic unit vector), ì  is the magnetic field,  = 0.5| ì |2 +  th +  CR + | ì | 2 /(8) is the total energy density,  th is the thermal energy density, and  CR is the CR energy density,  tot =  th +  CR +  mag is the 1 https://bitbucket.org/rteyssie/ramses/src/master/sum of the thermal  th = ( − 1) th , CR  CR = ( CR − 1) CR and magnetic  mag = | ì | 2 /(8) pressures, with adiabatic indexes  = 5/3 and  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 the anisotropic diffusion flux ( is the diffusion coefficient, assumed constant in this work), L rad = L rad,th + L rad,CR→th the total radiative losses which includes thermal and CR radiative losses, respectively 2 .The CR loss term encompasses both energy losses due to hadronic and Coulomb interactions L CR = −7.51× 10 −16 ( e /cm −3 )( CR /erg cm −3 ) erg s −1 (see e.g., Ensslin 2004;Guo & Oh 2008) transferred to the thermal component at a rate H CR→th = 2.63 × 10 −16 ( e /cm −3 )( CR /erg cm −3 ) erg s −1 .
In the CR energy equation (Equation ( 5)), ì  st = ( CR +  CR ) ì  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 selfexcited Alfvén waves in different gas phases (e.g.Armillotta et al. 2022), allowing for streaming velocities exceeding ì  A .As a preliminary step towards a deeper understanding of CR feedback in a cosmological context, we assume that CR streaming velocities cannot reach super-Alfvénic velocities, and leave the implementation of self-consistent CR diffusion coupling (like the one presented in Farber et al. 2018;Xu & Lazarian 2022;Armillotta et al. 2022) for future work.Furthermore, since CRs scatter on the Alfvén waves, they experience a drag force, transferring energy to the thermal gas at a rate The typical value of  ≃ 3 × 10 28 ( CR /1GV) 0.34 cm 2 s −1 depends on the particle rigidity  CR (Jóhannesson et al. 2019), which for 1 GeV proton has a value of  CR = 1 GV.Previous efforts to constrain  in CR simulations by comparing Fermi LAT -ray observations(e.g.Fire-2, arepo and ramses, Chan et al. 2019;Werhahn et al. 2021;Nuñez-Castiñeyra et al. 2022, respectively) have resulted in large discrepancies in the inferred diffusion coefficient favoured by observations.Nuñez-Castiñeyra et al. (2022) argue that, although differences in the assumed CR spectrum or -ray emissivity per atom can reduce the discrepancy by factors of up to 2.5, these are not sufficient to bring different simulations to agreement.They suggest that differences in the thermal state of the gas and the magnetic field properties, achieved by different MHD codes, could be the cause of these differences.The exploration of how different CR transport mechanisms and diffusion coefficients affect the -ray luminosity is the scope of an upcoming paper (Rodríguez Montero et al., in prep).We adopt here a value of 3×10 28 cm 2 s −1 favoured by the results of isolated CR simulations using similar versions of ramses to the one used in this work (Dashyan & Dubois 2020;Farcy et al. 2022;Nuñez-Castiñeyra et al. 2022), and leave the analysis of a higher  = 3 × 10 29 cm 2 s −1 (favoured by the Fire-2 simulations Chan et al. 2019;Hopkins et al. 2021b) to future work.

The nut suite of simulations
The new suite of simulations (see Table 1 for details) presented in this work uses the initial conditions of the nut suite (Powell et al. (2011) Figure 1.Overview of the CRMHD model of the nut simulation at  = 1.5.Background: Colour composite of the zoom-in region centred on the nut galaxy, with gas density (silver blue), gas temperature (orange) and CR energy density (yellow), showing 600 kpc along the x-axis.Three cold filaments are feeding the growth of the central galaxy, where strong feedback events are driving large-scale outflows of hot and CR-dominated gas.The effects of feedback are also observed in smaller objects embedded in the infalling filaments.Left inset: RGB colour composite zoom into the immediate galactic region as would be observed with the F200W, F150W and F090W filters from JWST.This mock observation is generated using the SKIRT radiation transfer code by assuming a fixed dust-to-metal ratio of 0.4 and the BARE-GR-S dust model from Zubko et al. (2004).Right inset: Edge-on view of the nut galaxy presenting different gas quantities modelled in the CRMHD simulation.In clockwise order from top-left corner: gas density, gas temperature, CR energy density, and magnetic field energy density.The cold and dense disk is shown in this projection, combined with a large hot outflow along the axis of rotation of the galaxy.CR and magnetic energies fill the halo and are shaped by a recent outflows.and Geen et al. (2013)) generated at  = 499 with cosmological parameters consistent with the results of WMAP5 (matter density Ω m = 0.258, dark energy density Ω Λ = 0.742, baryon density Ω b = 0.045, amplitude of the linear power spectrum at the 8 Mpc scale  8 = 0.8 and Hubble constant  0 = 72 km s −1 Mpc −1 Dunkley et al. 2009).The computational box has a side of comoving length 9 Mpc ℎ −1 with a spherical region of approximately 2.7 comoving Mpc ℎ −1 in diameter defining the volume in which high resolution is performed.Within this zoom-in region a  vir ( = 0) ≃ 5 × 10 11  ⊙ DM halo hosts a Milky Way-like galaxy (from this point simply identified as the nut galaxy) with a DM particle resolution of  DM = 5.5 × 10 4  ⊙ .To probe the multi-phase structure of the ISM, we allow refinement within this region to be triggered via a quasi-Lagrangian method, for which cells with total mass (i.e.gas, stars and dark matter) above 8 DM Ω m /Ω b are subdivided into 8 equal cells, reaching a maximum physical resolution of cell size Δ max ≃ 23 pc approximately constant in proper units (an extra level of refinement is triggered for expansion factors of  exp = 0.1, 0.2, 0.4, 0.8).
In Fig. 1 we show an overview of our fiducial simulation with CRs and MHD centred on nut.The main galaxy is fed by three large filaments (gas density in silver), and at their intersection there is a gas halo filled with high temperature gas (in orange) and CR energy (yellow) from recent feedback events.Additionally, we show zoomin insets of the galactic region of nut seen face-on (left circle) as it would be seen through JWST F200W (red), F150W (green) and F090W (blue) filters, and edge-on (right circle) density-weighted projections of the raw simulation data.The mock JWST RGB image has been obtained with the radiation transfer code skirt (Camps & Baes 2015), assuming a fixed dust-to-metal ratio of 0.4 and the BARE-GR-S dust model from Zubko et al. (2004).At  = 1.5 nut shows signatures of a large feedback event: hot bipolar outflows emerge perpendicular to the cold, dense and magnetised disk.These outflows clear denser and colder gas in the CGM, creating bubbles seen in the gas density, magnetic energy and CR energy.The face-on view through SDSS filters shows a late-type spiral with a nuclear cluster at its centre dominated by an old stellar population, and a disk filled with lanes and filaments of dense dust absorption and small regions of recent star formation.For all our runs, we include metal-dependent radiative cooling with interpolated cooling curves from cloudy (Ferland et al. 1998), down to ∼ 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  = 9 (Haardt & Madau 1996), with gas denser than  H = 0.01 cm −3 considered to be self-shielded from UV radiation with an exponential dependence on  H . Furthermore, despite the high resolution of our simulation, we lack the cooling modelling capabilities to resolve primordial (Population III) SN explosions in minihalos (e.g.Whalen et al. 2008;Wise et al. 2012) that would potentially raise the initial gas metallicity to ∼ 10 −3  ⊙ and hence increase cooling at high redshift.To avoid this issue, we initialise the simulation with a uniform metallicity floor of 10 −3  ⊙ .
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 magnetothermo-turbulent (MTT) star formation model presented in Kimm et al. (2017) and in Martin-Alvarez et al. (2020) for its MHD version.Star formation can only happen where the gas number density is above 10 cm −3 and the gravitational pull is larger than the combined support from turbulent, thermal, magnetic and CR pressure 3 .We determine the rate of gas to stellar mass conversion following a Schmidt law (Schmidt 1959) where we allow  ff to be a local parameter that is determined by the gas cell conditions as predicted by the multi-free-fall model from Hennebelle & Chabrier (2011) (generalised to the model of PN11 by Federrath & Klessen ( 2012)), and refer the reader interested in further details of this model to appendix B in Martin-Alvarez et al. (2020).
Ultimately, the conversion of gas to stars is modelled by sampling stochastically a Poisson mass-probability distribution for the final stellar particle mass (Rasera & Teyssier 2006), with the minimum stellar mass at creation fixed at 4.5 × 10 3  ⊙ .
As discussed in the introduction, instead of performing a calibration process to determine artificial 'boosting' factors to the canonical type II SN energy of 10 51 erg that may allow, e.g., the recovery of a stellar-to-halo-mass (SMHM) relation in agreement with local observations (Rosdahl et al. 2018), we model SN feedback using the 'mechanical' approach of Kimm & Cen (2014).Considering the local gas conditions and the available spatial resolution, it determines whether the Sedov-Taylor expansion phase of the SN can be captured directly (Blondin et al. 1998;Thornton et al. 1998) or if numerical overcooling will artificially affect its evolution.In the latter case, momentum prescribed by the snowplow phase is injected.In order to determine when a stellar particle should experience a SN feedback event, we follow the multiple explosions method by Kimm et al. (2015), in which we use the real delay of SNe given by the population synthesis code starburst99 (Leitherer et al. 1999), with SNe taking place as early as 3 Myr and as late as 50 Myr.These SN events inject mass, momentum and energy into the host and a maximum number of 48 neighbouring cells, with a canonical efficiency of  SN = 10 51 erg/(10 ⊙ ) 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 ì  field.Several mechanisms have been presented over the last decades as candidates to generate galactic magnetic fields of the order of 10 −7 − 10 −5 G (e.g.Mulcahy et al. 2014;Beck 2016).These can be classified by their cosmological (e.g., electroweak or quantum chromo-dynamics phase transitions Vachaspati 1991;Olesen 1993) or astrophysical (plasma seeding processes like Biermann battery Biermann 1950;Pudritz et al. 1989) nature (see Durrer & Neronov 2013;Vachaspati 2021, for recent reviews on the matter).Whether magnetic fields in galaxies are primordial in nature or rapidly amplified by astrophysical or dynamo processes (Pakmor & Springel 2013;Rieder & Teyssier 2017;Vazza et al. 2017;Martin-Alvarez et al. 2018, 2021), the theoretical expectation is that they reach ∼ 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  axis4 .In order to mimic such rapid amplification, we employ a uniform primordial magnetic field of strength  0 = 3 × 10 −12 G.This simple magnetic field initialisation directly provides the expected strengths of 1 − 10  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  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 5 a new cosmological simulation analysis package which allows for a smooth interfacing of halo and galaxy selection codes, in-depth analysis algorithms and simulated mock observations.It is developed with an easy-access philosophy in mind, combining Python for accessing the data catalogues and Fortran 90, openMP parallelised routines to access the raw ramses data.It takes a cosmological simulation snapshot and outputs a convenient galaxy and halo catalogue following the HDF5 6 format standard, which makes possible the study of many physical properties in the simulation without the need to read back again the original binary files.Additionally, these files are usually less than 0.3% the size of the original outputs, making simple analysis and comparison between simulations much less memory and computationally expensive.Ozymandias has been written to work with two different halo and galaxy finders (although the I/O for halo catalogues is fully compatible with other halo finders, providing the data structure information): the HaloMaker code by Tweed et al. (2009), a 3D density-threshold structure finder extended to use the shrinking spheres centring method (Power et al. 2003); and the VELOCIraptor code (Elahi et al. 2019;Cañas et al. 2019) in its enhanced version by Rhee et al. (2022), a robust, scale-free, phase-space structure-finder.Once halos and galaxies are identified, a hierarchical structure is built between them, linking galaxies with their host halos and linking satellites with central galaxies.Then, in order to allow for an easy tracking of each galaxy/halo across cosmic time, particle lists for each identified object are compared between consecutive snapshots, allowing for an easy access to the progenitors Table 1.Runs of our new nut suite of simulations studying CR feedback included in this work.Columns indicate from left to right the simulation name, their employed solver and the details of the CR physics if CRs are present.All simulations use the same initial conditions, a MTT star formation prescription and the canonical configuration for mechanical feedback prescription.For simulations using the MHD or the CRMHD solver, all simulations employ a uniform comoving initial magnetic field of  0 = 3 × 10 −12 G.

Simulation Solver CR physics HD
CRMHD Streaming and  = 3 × 10 28 cm 2 s −1 in the previous snapshots.We use this package to identify the main galaxy, and refer to the galactic region as the sphere with radius 0.2 vir,DM ,  vir,DM being the host DM halo virial radius.Quantities such as the galaxy stellar mass  * and gas mass  gas are measured within this spherical volume at every snapshot.

The global effect of cosmic rays in the evolution of nut
In Fig. 2 we show density-weighted projections of nut at  = 1.5 for the three simulations (HD, MHD and CRMHD), comparing faceon 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  = 1.5, the nut galaxy has formed a thin gas disk in all simulations, that has large-scale spiral features with high density clumps and low density inter-arm regions.In the HD and MHD simulations (top and middle rows) the galaxy has the general properties of a disk galaxy within a virialised halo: a cold ( < 10 4 K) thin disk with close to solar metallicity Z ⊙ is embedded within a low density, hot ( > 10 6 K) CGM.The edge-on views aid in identifying long filaments of cold, dense and low metallicity gas inflowing into the galaxy.In the CRMHD simulation (bottom rows), while the dense and cold region of the disk can still be clearly seen in the projections, there are significant differences with respect to the no-CR simulations.Firstly, the gas disk shows much more filamentary substructure, and in comparison to the HD and MHD simulations the maximum value of density in the clumps is substantially reduced, as well as the volume occupied by low density, hot gas in the inter-arm regions.
Secondly, the high density centre of the galaxy in HD and MHD extends to ∼ 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  ∼ 10 5 K. Against this lower temperature 'background', recent hot outflows can clearly be distinguished in the edge-on view, although the lower average metallicity of the disk gas compared to HD and MHD makes the metal enrichment of outflows less prominent.Instead, the CGM of nut is now more uniformly enriched to ∼ 0.1 Z ⊙ .Furthermore, the mock JWST images of the HD and MHD simulations show nut as a disk galaxy dominated by a massive bulge and an old stellar population, with dark, dense filaments seen as dust lanes in the edge-on views.Also, massive stellar clumps formed by largescale disk instabilities persist due to inefficient feedback (e.g.Inoue et al. 2016;Bournaud et al. 2013).In comparison, the CRMHD simulation shows a dramatic bulge reduction, with the central luminosity of the galaxy reduced by 0.5 dex.Additionally, instead of having massive clumps within its disk, small star forming regions are seen all across the disk, with the luminosity of young stellar population more extended.To better quantify this important difference, we have measured the compactness of the stellar emission of nut using the total luminosity in the JWST filters.By fitting a 2D Sérsic profile, we obtain half-light radii  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): with a larger value of log Σ 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  * and  50 are comparable.When considering the CRMHD simulation, a significantly reduced  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 panel) and gas (right panel) masses inside the galactic region are shown for the HD model with blue lines in Fig. 3.The stellar mass growth follows the common trend for a MW-like halo (e.g.Kim & Ostriker 2015;Martin-Alvarez et al. 2018): 1) accretion phase: significant gas accretion at high redshifts allows for a fast growth of stellar mass, reaching 10 10 M ⊙ at ∼ 1.5 Gyr after the start of the simulation; 2) feedback phase: feedback from SNe and the formation of a virial shock after  ∼ 4 slows down star formation, leading to self-regulation with  * ∼ 4 × 10 10 M ⊙ at  = 1.5, in agreement with the equivalent MFBm model in Kimm et al. (2015).The evolution of  gas is also in agreement with the findings from Kimm et al. (2015), exhibiting a rapid growth until  ∼ 5. Subsequently, the fast depletion timescale reduces the gas mass inside the galactic region until  ∼ 3. A new inflow of gas mass takes place from  = 3 to  = 2, coinciding with the last major merger that the nut galaxy experiences (Martin-Alvarez et al. 2018, e.g.).Fig. 3 includes the MHD run with a primordial magnetic field of 3 × 10 −12 G but no cosmic rays.This illustrates that the moderate magnetisation employed here does not significantly affect the global properties of the galaxy (Su et al. 2017;Pillepich et al. 2018;Hopkins et al. 2020;Martin-Alvarez et al. 2018).
However, the addition of CRs gives rise to a radical difference in the star formation history (SFH) of our nutgalaxy: the steepness of the stellar mass growth is reduced during the accretion phase, resulting in  * ≃ 2 × 10 9 M ⊙ , ∼ 1 order of magnitude lower with respect to the HD and MHD models at  ∼ 4.During this period of high-redshift gas accretion, the CRMHD run has only a ∼ ×2 − 4 deficit in  gas inside the galactic region.Therefore the decrease in  * cannot be solely attributed to the depletion of  gas from the galactic region.We also include on the right panel the cold gas mass (dashed lines) as defined by the constant entropy limit7 of Gent (2012), showing that for the CRMHD simulation the cold gas mass decreases by ∼ 1 dex compared to HD and MHD at  ≳ 4.This points at the depletion of the cold gas mass in the ISM as the major cause of the lower stellar mass growth at these cosmic times.Following this period of strong star formation suppression, the CRMHD run exhibits a similar period of gas accretion beginning at  ∼ 3 as the 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  50 obtained by fitting the RGB images with a 2D Sérsic profile.For each simulation the stellar mass within 0.2  vir,DM and galaxy compactness log Σ 1.5 (Barro et al. 2013), is listed, with a larger value indicating a more compact stellar emission.Compared to HD and MHD simulations, the disk of nut shows much more structure in the CRMHD simulation, dominated by a more diffuse ISM with less compact and dense clumps.The CGM is also strongly altered, with the cold dense filaments of HD and MHD replaced by a diffuse and cold halo more uniformly enriched by metals.The mock observations of HD and MHD show a galaxy dominated by an old stellar population and a large bulge, while the nut galaxy in CRMHD is transformed into a less compact system dominated by a disk.In this case, star forming regions are scattered across spiral arms and clumps all over the disk.The right panel distinguishes the total gas mass (solid line) from the cold gas mass (dashed line), as defined by the constant entropy limit of the cold phase of the ISM by Gent (2012).In all runs stellar mass grows rapidly up to  ∼ 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 ∼ 1 dex at  ≳ 4. Below  ∼ 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  = 1.5 to be similar to that in the HD and MHD runs, whilst the stellar mass remains 4 times lower.rest of the runs, which eventually translates into a new phase of steep  * 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 ∼ 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 ( * ) to halo mass ( halo ) relation for the HD (blue diamonds), MHD (purple stars) and CRMHD (green triangles) runs at four different redshift ranges from  = 1.5 to  = 8.We define the halo mass as the full mass enclosed in gas and DM within  vir,DM and the stellar mass of the galaxy.These redshift ranges are separated by vertical dashed lines with their central redshift values indicated at the bottom of the panel.Results from the abundance matching models by Behroozi et al. (2013) (orange shaded region) and Moster et al. (2013) (blue shaded region) at  = 1.5 are included for reference.The HD and MHD runs show very similar behaviour across the redshift range presented, with only a small stellar mass suppression by magnetic fields.Both feature large baryonic conversion efficiencies, close to the line of total conversion of gas into stars for a given halo, given by  b  halo , where  b is the universal baryon fraction defined as Ω b /Ω m , the ratio of the cosmological baryon and matter density parameters.This result is consistent with the well established understanding of non-calibrated, canonical SN feedback: it fails to reproduce the results of abundance matching models, leading to a too massive stellar mass for a given dark matter halo mass.The HD run begins with ∼ 32% baryon conversion efficiency and grows to ∼ 56% by  = 2, while the MHD run has a slightly lower efficiency of ∼ 47%.The CRMHD run shows a consistent decrease of the baryonic conversion efficiency across all redshifts, with a maximum value of ∼ 14% at  = 2 and a minimum of 5% at  = 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  = 1.5.Despite the decrease of almost an order of magnitude in  * for the CRMHD run compared to the HD and MHD runs, the CRMHD model does not follow the 'natural' extrapolation of the abundance matching models to lower halo masses.It is important to point out that reliable constraints on the number density of galaxies hosted by low mass halos of 10 10 M ⊙ between  = 6 and  = 8 are scarce, adding uncertainty in abundance matching models at high redshift.Additionally, the observational stellar mass functions used in these works are subject to observational incompleteness below 10 7 M ⊙ (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.

How cosmic rays alter the gas distribution and the 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  = 8 and  = 4.This epoch also is where we find the largest discrepancy in  * between the CRMHD and HD/MHD runs.However, the 10 times smaller stellar mass measured at these redshifts cannot be solely attributed to the lower gas content in the galaxy which is much more moderate, and we need to consider whether the star formation efficiency of the gas is reduced compared with the HD and MHD runs.To explore the physical conditions of the  ISM gas in nut, Fig. 5 shows temperature-density phase diagrams for the three runs used in this study.Phase diagrams based on a single simulation snapshot can be hard to interpret as they are subject to the stochasticity driven by processes such as feedback events and mergers.To avoid this issue, the phase diagrams in Fig. 5 are averaged over multiple snapshots within the galaxy's dynamical time at  ∼ 2.5, weighted by their difference in time to  ∼ 2.5.We estimate the galaxy's dynamical time  dyn as the time required for a test particle to complete one full orbit at  = 0.2 vir,DM with circular velocity: where  (< ) =  b +  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  dyn of the three runs compared, which for  ∼ 2.5 corresponds to the CRMHD simulation and has a value of  = 0.556 Gyr.Therefore, our measurements have a central redshift ⟨⟩ and a physical time range .
For the HD run, the left panel of Fig. 5 shows that the low density gas (10 −26 g cm −3 to 10 −24 g cm −3 ) in the galactic region is dominated by gas with temperature  ∼ 10 4 K (where hydrogen recombination and resonance lines of atomic hydrogen dominate the cooling of gas) and above, a typical signature of SNe shock-heating.Gas accumulates in the  ∼ 10 4 K region due to thermal instability (Field et al. 1969;Wolfire et al. 1995).When this gas is perturbed by processes such as shock compression, it collapses into cold dense clumps that populate the tip of the phase diagram in the bottom right corner.In agreement with the results found by Martin-Alvarez et al. (2018, 2020) for nut, the addition of magnetic fields helps sustain a somewhat larger fraction of gas in hotter phases, as well as allows for gas cells to reach lower densities whilst retaining their temperature.This is reflected in the size of the 'dark' regions of the phase diagram with densities below 10 −24 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 : cold and dense (blue region), warm (orange region) and hot and diffuse (red region).This separation is inspired by the work of Gent (2012).Similar amounts of gas are found in the hot and warm regions of the HD and MHD runs, with the MHD run having a slightly higher gas mass in the warm phase (from 49% in HD to 55% in MHD).The hot phase only contributes 3% to the total gas mass in both HD and MHD.
However, the gas phase diagram for the CRMHD run (right panel of Fig. 5) shows marked differences compared with the non-CR models: (a) at densities  < 10 −24 g cm −3 , gas is no longer confined to temperatures above 10 4 K when  < 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. < 10 4 K with 10 −24 <  < 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 two compared to the HD model as this gas now resides in the warm medium.Notably, it is in these dense regions of the phase diagram that we expect gas to undergo a more efficient star formation.This considerable reduction of cold gas mass explains, at least partially, the causal link with the measured decrease of stellar mass in this simulation as inferred from Fig. 3.Note that we have not yet quantified the turbulent or magnetic support of the gas, and how these affect the efficiency of star formation.Indeed, in combination with the thermal pressure, the additional CR pressure is able to support gas against further collapse, and the build up of a pressure gradient allows for the warm gas reservoir above the disk to be kept there for longer (e.g.Dashyan & Dubois 2020;Chan et al. 2022;Farcy et al. 2022;Armillotta et al. 2022).This depletion of cold gas, now transferred to the warm phase in the presence of CRs is in good agreement with other works (e.g.Dubois et al. 2019;Nuñez-Castiñeyra et al. 2022).By considering the medians of the 2D phase diagrams, two main features are readily visible: (i) the larger contribution to the mass budget by gas between 10 −24 g cm −3 to 10 −22 g cm −3 in the CRMHD run creates a more pronounced concave inflection point towards lower  at  ∼ 10 −23 g cm −3 ; (ii) gas mass 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  ≥ 10 5 K is greatly reduced, particularly at higher densities.

Distribution of supernova sites and feedback effectiveness
The modifications we observed to the gas distribution across ISM phases also have an effect on the distribution of SN explosions across the multi-phase ISM.These local properties will affect the efficiency of SN feedback on self-regulating star formation and modifying the ISM properties.We explore this issue by recording the local ISM conditions for each SN event.In Fig. 6 we plot the distribution of SN site gas density for two redshift ranges: between 3 <  < 6 (top panel), i.e. the period of largest  * difference between CRMHD and HD/MHD, and between 2 <  < 3 (bottom panel), which samples the period of fast growth of  * in CRMHD.We find for the HD simulation that (as in previous works Kimm et al. 2015) these PDFs are intrinsically bimodal with the two peaks at low and high densities, characterising the hot and cold phases of the ISM (see e.g.Gent 2012).However, the high density peak dominates the overall distribution, which at higher redshifts (top panel) resembles more a log-normal distribution of the cold ISM.When SNe explode in the HD case, young star particles have not had the chance to leave the dense regions of the galaxy in which they were born.This means the 10 51 erg from a single SN 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 (∼ 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 super-bubbles.
The presence of this bimodality in SN sites is strongly altered for the CRMHD simulation compared to no-CR runs.Although the 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 Sedov-Taylor phase.
effects seen in the MHD simulation, which arise from the support of magnetic fields, should also be in place in the CRMHD simulation, we have already seen (in Fig. 5) that the presence of CR pressure and heating mechanisms dominates over the effect of magnetic fields.The ISM in the CRMHD run is strongly depleted of high density cold gas, with the high density peak of the bimodal distribution now at ∼ 10 −23 g cm −3 in the CRMHD run instead of ∼ 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 (≤ 10 −22 g cm −3 ), with a large number of SNe exploding at  ∼ 10 −27 g cm −3 .We attribute this to star formation happening in a less clustered form throughout the galactic region, with the consequent feedback events being able to sample lower density regions shaped by CR pressure and previous SN-driven bubbles.In addition, SN remnants evolve for longer and maintain higher temperatures past the canonical evolutionary sequence (e.g.Naab & Ostriker 2017) when the effects of CR pressure are taken into account (see Rodríguez Montero et al. 2022 for a parameter space study of CR feedback effects on the evolution of individual SN using 3D MHD simulations).Furthermore, similarly to the behaviour of fully thermal SN remnants (Thornton et al. 1998), when CRs are accounted for in the evolution of SN explosions in low density environments, they can reach a larger value of deposited momentum compared to high density environments (Rodríguez Montero et al. 2022).
The bottom panel of Fig. 6 shows the PDF of SN site densities instead between the redshifts of 3 and 2. In this lower redshift range we find that the CRMHD PDF begins to converge to the ones obtained for the HD and MHD runs, although maintaining a higher fraction of SNe happening at lower density (i.e.∼ 10 −27 g cm −3 ).nut transitions into a rotationally supported (i.e.ratio of rotational velocity to dispersion velocity larger than unity) disk at  ∼ 3.5 (Martin-Alvarez et al. 2020) in HD and MHD simulations, while we found that for our CRMHD run this transition occurs closer to  ∼ 2.3.The nut galaxy 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 ∼ 10 −27 g cm −3 and a peak at ∼ 10 −23 g cm −3 .SN events taking place at lower densities during this period are more effective in suppressing star formation.This bimodality for CRMHD is still present in the lower redshift range (bottom panel), but the contribution of the low density peak is significantly suppressed, with the overall distribution resembling more the one obtained for the HD and MHD simulations.
is dominated by the gas velocity dispersion between  = 6 − 3, which means that CRs released by SN events within the galactic region can easily diffuse across the ISM and affect additional star-forming regions.However, in the lower redshift range of the bottom panel we are instead capturing the morphological transformation into a dynamically cold disk.With a disk morphology in place, CRs will be more prone to escape perpendicular to the disk, having a less drastic effect in the cold phase of the ISM at lower redshifts.We predict a dual behaviour depending on the morphology of the galaxy, which we will explore in future work.

Influence of cosmic rays in the launching of outflows
The CRMHD simulation has both lower  gas and  * at high redshift than the HD and MHD models.This must be, at least partially, due to differences in the outflow and inflow rates of these simulations.
In this work we focus on the analysis of outflowing gas.We measure instantaneous outflow rates  outflow as where we sum over all cells with positive radial velocities   , in a thin spherical shell of width  shell at a distance  from the galactic centre.The width of the shell is taken to be 0.01  vir,DM , which provides a consistent measure across simulations and cosmic time.Each cell is characterised by its gas density   , the radial velocity   , and the cell width Δ  .Previous studies with the nut galaxy (Tillson et al. 2015) found the contribution to the inflowing gas from satellites to be small (∼ 1%) compared to that from filaments.Nevertheless, for the analysis of the outflowing gas, satellites that cross our thin shell may contain gas that could appear to be flowing radially outward in the frame of the central galaxy (either due to satellite rotation, or because a given satellite crosses the measuring shell on its way back to the outer halo regions after its apocentric passage).To avoid these subtleties, we subtract all contributions in the thin shell from passing substructures.We thus remove gas located within the tidal radius of DM sub-halos containing more than 1000 particles.The tidal radius is estimated as the Jacobi radius with an isothermal sphere approximation for DM halos given by Eq. 8.107 of Binney & Tremaine (2008).Finally, we divide outflowing gas into separate temperature phases as follows: • Hot gas:  > 10 5 K, • Warm ionised gas: 10 5 >  > 9000 K, • Warm neutral gas: 9000 >  > 1000 K and • Cold gas:  < 1000 K .
In Fig. 7 we show the results of this outflow analysis in a thin spherical shell at 0.2  vir,DM for the three simulations across ∼ 4.5 Gyr of evolution, where  outflow (top panel), star formation rate (SFR) averaged over 100 Myr (middle panel), and the mass loading factor  =  outflow /SFR (bottom panel) is shown.Solid lines are the result of a running median smoothing with a size of 200 Myr.Note that we do not restrict our analysis to particular assumptions about the possible time delay between a star formation event in the galaxy and the consequent outflow at 0.2  vir,DM .We distinguish three distinct time phases when comparing the HD, MHD and CRMHD simulations: i) between  ∼ 13 and  ∼ 7 CRMHD reaches an outflow rate of 2 − 4 M ⊙ yr −1 , whereas the HD and MHD runs are closer to 0.1 − 0.4 M ⊙ yr −1 ; ii) between  ∼ 7 and  ∼ 3 the outflow rates of the HD and MHD runs reach 1 M ⊙ yr −1 (with the MHD run having overall larger values than HD), during times of strong starbursts.However, given their average SFR of 10 M ⊙ yr −1 , their mass loading factors remain well below unity, while the CRMHD simulation consistently exceeds  ∼ 1 during starbursts; and iii) between  ∼ 3 and  ∼ 1.5 the SFR of the three runs peaks at  ∼ 2.5 at a value of 10 M ⊙ yr −1 (and slowly drops to 2 − 4 M ⊙ yr −1 ).After this time, only CRMHD maintains a steady outflow rate of 8−9 M ⊙ yr −1 while MHD drops to 2 M ⊙ yr −1 and HD is kept at 2 − 4 M ⊙ yr −1 , with a jump to 8 M ⊙ yr −1 around  ∼ 2. It is also important to note that in comparison with the HD and MHD runs, the CRMHD runs experiences a fast rise in SFR of 1 order of magnitude between the redshifts of 3 and 2.3.As mentioned before, during this period nut experiences the merger with a large companion, causing a large increase in the cold gas mass of the galaxy (see Fig. 3) and a burst in star formation.We find that this period is connected with a large increase of inflow rate to nut in the CRMHD run (Rodríguez Montero et al., in prep.), caused by the cold, high density extraplanar gas supported by CR pressure to become unstable during the merger interaction.Overall, the cosmic evolution of  outflow appears more steady than the SFR for the CRMHD run.This suggests that while SN feedback is rarely efficient in the HD and MHD runs with large amount of gas consumed into stars, in the CRMHD run large quantities of gas are ejected from the galaxy and gas-to-stars consumption rate is much lower.
In order to further understand the underlying differences between outflows in the non-CR runs and CRMHD, we compare their outflow Bottom panel: Mass-loading factor  computed using the running median results of  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 all cosmic time explored.Given the large suppression of star formation for  ≥ 3, the CRMHD simulation has outflow mass-loading consistently exceeding unity, only reached by the HD and MHD simulations during the merger-induced starburst at  ∼ 3 − 2. velocity distributions during two distinct phases of the SFH: 1) a peak in SFR present in all simulations at  ∼ 4; and 2) a period of similar SFRs across all runs (∼ 8 − 10 M ⊙ yr −1 ) occurring approximately at  ∼ 2.5.The latter is characterised by a merger with a small companion.We compute stacked mass-flow rate histograms following the same method detailed in Section 3.2 for weighted stacking of PDFs, employing snapshots within a dynamical time of the galaxy.These histograms show the total  separated into radial velocity bins for the same thin shell employed in Fig. 7 at 0.2  vir,DM .We show the HD, MHD and CRMHD simulations in blue, purple and green, respectively.The top panel of Fig. 8 shows that for  = ⟨3.94⟩ the outflowing gas carrying the majority of the mass outside of the galaxy is found at lower radial velocities in the CRMHD run than in HD and MHD.While the HD distribution is dominated by gas between and CRMHD (green) simulations, stacking the histogram values in snapshots spanning a  dyn range, using the method detailed in Section 3.2.Vertical dotted lines show the stack-weighted averages  esc for the redshift quoted, with the corresponding standard deviations denoted by the shaded regions.Values for the total mass outflow rate above  esc (  esc ) are included in each panel for each simulation.The peak of the outflowing gas distribution is both more prominent and displaced to lower radial velocities in the CRMHD simulation in both redshift ranges compared to HD and MHD.At high redshift, the shallower gravitational potential in the CRMHD simulation allows a much larger fraction of the high velocity tail to exceed the local  esc .
60 km s −1 and 120 km s −1 , it rapidly drops towards lower velocities, with a small fraction of the total outflow rate (0.05 M ⊙ yr −1 ) exceeding the  esc of the galaxy (given by the vertical dotted blue line).This escape velocity is computed for each time  by counting the total mass enclosed within 0.2  vir,DM .The MHD simulation instead shows a tail that extends to higher velocities, contributing 10 times more to the mass outflow rate above  esc than HD.However, despite this higher efficiency driving fast outflows, the majority of the outflowing mass between 60 km s −1 and 120 km s −1 remains the same, and the mass above  esc is 3 times lower than the CRMHD.This is due to the fact that  esc is half the value for CRMHD than for HD/MHD, a signature of a shallower potential well.This is due to the CR feedback significantly affecting the growth of the central stellar mass of nut at high redshift.Therefore, the larger outflow rates of CRMHD at  ∼ 4 compared to HD and MHD can be explained by a larger fraction (almost 1 dex compared to HD and MHD) of 'low-velocity' gas and an easier escape of fast outflows in a less deep gravitational potential.To further explore these trends, we also analyse the outflow velocity distribution at  ∼ 2.5 (see bottom panel of Fig. 8).Although now the fast growth of nut in the CRMHD simulation between  = 3 and  = 2 somewhat bridges the gap in  esc , this simulation still has 6 times as much gas outflowing above  esc than the HD and MHD runs.This is not surprising given that during this period  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  vir,DM lacks the kinetic energy necessary to escape the local gravitational potential in all simulations.However, it remains an open question whether this gas needs to have the energy necessary to escape the halo in order to reduce the reservoir of gas available for star formation or if it is enough to prevent its re-accretion over sufficiently long timescales (see e.g.Somerville & Davé 2015, for preventive feedback in the context of AGNs).As far as our current analysis goes, the CRMHD model is able to reach higher mass loadings than the no-CR runs.These outflows are mainly comprised of a slow moving wind that, without any further acceleration beyond 0.2 vir,DM , will not be able to escape the halo and will eventually fall back into the galaxy if undisturbed by future feedback events.

Dynamics of multi-phase 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 disk.To provide evidence for this higher efficiency causing these higher mass-loading results in CRMHD, Fig. 9 shows the  outflow and  outflow -weighted radial velocities of outflowing gas  outflow in the HD (blue), MHD (purple) and CRMHD (green) simulations, separated in the different phases we explore.HD and MHD outflows are dominated by the hot and ionised phases during the majority of the evolution of nut.
Outflowing in the warm neutral and cold phases in these runs are identified with sudden, low velocity peaks in the velocity plots, which we found to be related to clumps of gas in eccentric orbits around nut, rather than genuine outflows.Instead, the CRMHD run shows that the hot phase never dominates the outflowing mass rate.The warm ionised phase dominates (1 − 10 M ⊙ yr −1 ) the outflowing gas mass up until  ∼ 4, with considerable contributions of warm neutral and cold gas at higher redshifts.We find that multi-phase outflows (∼1:2:4:2 ratios for cold, warm neutral, warm ionised and hot phases, respectively) are only found consistently in the CRMHD simulation, and at high redshift ( ≥ 4).The HD and MHD runs fail to produce such gas mass distributions across phases at all redshifts presented here.We compare the radial velocities with the escape velocity at 0.2  vir,DM .To understand whether the outflowing gas is gravitationally bound, we compare these velocities with the local escape velocity shown in the four panels with dashed lines.While a radial velocity higher than  esc is not a sufficient condition for gas to completely escape the gravitational potential of a halo, it is a useful proxy commonly used in observations and theoretical studies.The radial velocities of the hot phase (top left panel) range from ∼ 100 km s −1 for  ≥ 3 to ∼ 80 km s −1 at lower redshifts in the no-CR simulations.For the case of CRMHD, the velocities are consistently higher, being in the range ∼ 100 − 200 km s −1 .These outflow-weighted velocities are in agreement with observations of the warm ionised gas in local dwarf galaxies by McQuinn et al. (2019) and Marasco et al. (2023), while falling short of the ∼ 400 − 600 km s −1 value reported in Concas et al. (2022) for their sample of galaxies between 1.2 <  < 2.6.However, the CRMHD simulation is the only one consistently able to reach outflow velocities above the escape velocity within the galactic region, up until  ∼ 2.5, when its distribution becomes similar to the non-CR cases.This behaviour is in agreement with the histograms shown in Fig. 8.
Outflows in the warm neutral and cold phases are rarely found in the HD and MHD runs (except during the merger taking place at  ∼ 3 − 2), and only appear as rare and sharp peaks in the   plots.By examining projections of the galactic region of nut at the redshifts corresponding to these peaks, we have found that for HD and MHD all cold and warm neutral outflows are linked to extended tidal features of satellites, and are not spatially related to the rest of hot and warm ionised outflowing gas.In fact, CRMHD is the only run that shows warm neutral and cold outflows across the majority of snapshots.These outflows are also spatially correlated with the hotter phases of the outflowing gas and are not dynamically linked with tidally stripped satellites (see Section 3.3.4and Fig. 11 for further details).Contrarily to the hot and warm ionised phases, the warm neutral and cold gas appear to show an evolution in outflow velocity, increasing from ∼ 50 km s −1 at high redshift to a maximum of ∼ 100 km s −1 at  ∼ 2.5 − 2. Similarly to the warm ionised phase, the median   of these two phases is below  esc (except for peaks at ∼ 2.2 and ∼ 2.5 Gyr in the cold phase).Overall, we find that these colder outflows struggle to reach the escape velocities at 0.2  vir,DM , even in the CRMHD run and despite its shallower gravitational potential.

Cosmic ray-driven re-acceleration of outflows in the CGM
The CRMHD simulation displays considerably higher outflow rates at 0.2  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  vir,DM in Fig. 10.This is carried out by computing the mass-flow-rate-weighted median of the radial pressure support for the outflowing gas in the thin shell at 0.2  vir,DM .As per the hydrostatic equilibrium equation, the negative pressure gradient of the gas needs to balance the inward pull of gravity for a parcel of gas not to inflow.The gravitationally-dominated regime of which this equilibrium situation is a limiting case, is indicated by the grey band below a value of 1.For a pressured fluid , this corresponds to − r • ì ∇  = r • ( ì ), where r is the radial unit vector, ì  is the local gravitational acceleration and  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,  ther , as solid circle symbols coloured in red linked by dashed lines, and CR pressure,  CR , as green circles connected by solid lines.
We begin by comparing the support for the hot phase, where we separate for clarity the HD (top left), MHD (middle left) and CRMHD (bottom left) simulations.We found that the hot and warm ionised phases show a large variance for these quantities, so we include the second and fourth quartiles of the underlying distribution as a shaded band in order to exemplify this.The three runs display a large amount of thermal pressure support against gravity for  ≥ 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.for the HD (blue), MHD (purple) and CRMHD (green) runs measured at 0.2  vir,DM .Each pair of panels show the instantaneous outflow rate (top) and the  outflow -weighted radial velocity of outflowing gas (bottom).The  outflow panels show the same quantities as Fig. 7: low opacity histogram is the instantaneous value, while the solid lines show the running median with a filter of 200 Myr.In the   panels, the escape velocity of each model at the shell radius is showed with dashed lines, for ease of comparison to the median  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 −1 for  ≥ 3 and ∼ 70 km s −1 for  ≲ 3, dominating the outflow rate in HD and MHD simulations.The hot phase in CRMHD reaches the largest outflow velocities (∼ 150 − 200 km s −1 for  ≥ 2.5 and ∼ 100 km s −1 for  ≲ 2.5), with a non-negligible fraction of the outflowing gas mass exceeding  esc consistently up to  ≥ 2.5.The dominant mass-flow rate is coming from the warm ionised phase in the CRMHD, while above  ∼ 5, warm neutral and cold phases largely contribute to the total outflow rate.Warm neutral and cold phases are scarce in the HD and MHD simulations, usually identified as low velocity peaks which can matched with crossing gas features of orbiting satellites.However, in the CRMHD run these outflows while slower than the hot gas can still reach ∼ 100 km s −1 at lower redshifts.
This change in the support of hot outflowing gas in all simulations is attributed to the growth of the virial shock at lower redshifts, as the cooling rate becomes smaller than the post-shock compression rate (Dekel & Birnboim 2006) and cold accretion becomes less important (Powell et al. 2011).A hotter CGM means a larger thermal pressure confining gas to the galactic disk, and buoyancy forces present at high redshift becomes less important.Although, this behaviour is observed for all runs, the CRMHD run features a larger frequency of radial support peaks exceeding the force balance threshold, despite the higher density of the CGM in the CRMHD simulation compared with the non-CR simulations (see Section 3.1).We attribute this to a colder CGM with a lower fractional contribution from thermal support due to the additional support of CR pressure (see Section 3.1).We do not include the contribution of magnetic support, as we find it to be negligible for outflows at all times for the MHD and CRMHD simulations9 .Thermal pressure dominates the hot phase in CRMHD .Pressure gradient support against gravity for outflowing gas separated into different temperature phases for HD (blue), MHD (purple) and CRMHD (green) simulations.The dividing line between the shaded and unshaded region indicates marginal support.In the shaded region, gravity dominates the dynamics, while in the unshaded region, gas is accelerated outwards by the pressure gradients.Outflowing gas is measured in a thin shell at 0.2  vir,DM and divided into the phases described in Section 3.3.2according to temperature.Shaded, coloured bands represent the second and fourth quartiles of the underlying distribution.The left column corresponds to the hot gas, with each row corresponding to the HD (top), MHD (middle), and CRMHD (bottom) simulations, respectively.While below  ∼ 3 the thermal pressure struggles to accelerate gas in the HD and MHD simulations, CR pressure consistently supports and accelerates the hot gas across all cosmic time.Only CR pressure can mildly support the warm ionised phase between 3 <  < 6, but it supports and accelerates warm neutral and cold gas for  ≤ 5. down to  ∼ 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  < 3.This change in behaviour coincides with the drop in hot gas radial velocity after  ∼ 3 shown for CRMHD in Fig. 9.While hot fast outflows in the HD and MHD simulations find it harder to travel through a the CGM at lower redshifts, in the CRMHD the additional support of CRs injected by SNe to the wind bubble further accelerates them when their thermal buoyancy stalls.
We compare in the right panels of Fig. 10 the pressure support against gravity in the warm ionised, warm neutral and cold phases, field injection by SN remnants (e.g.Martin-Alvarez et al. 2021), which we do not explore in this study.now for the three simulations combined.As warm and cold outflows are scarce in the HD and MHD runs, we limit our analysis for these phases to the CRMHD run.These confirm that the thermal pressure support (as well as the magnetic, which we do not include) in the HD and MHD runs is negligible.This further supports the idea that for the runs without CRs, outflowing gas colder than 10 5 K does not experience further acceleration at 0.2  vir,DM .Due to outflow velocities lower than  esc , outflows likely rain back down onto the galaxy as a galactic fountain.The CRMHD run features warm ionised outflows that are marginally supported against gravity, following an overall trend similar to that of the hot gas, but now with CR pressure dominating over thermal pressure at all times.In fact, CR pressure support dominates over thermal pressure for the three phases in the right column of Fig. 10.The warm neutral and the cold components of the outflow show large accelerations due to CR pressure, exceeding by up to 7 times the gravitational acceleration.These phases display a large variability, also found in the analysis of  r (see Fig. 9).This can be attributed to a combination of outflow 'burstiness' and a lower volume-filling of these phases.Therefore, whenever outflows entraining warm neutral and cold gas are present, they are usually subject to large CR pressure gradients that not only balance gravity but provide additional acceleration.This means that, although the outflow velocities of these two phases at these radii are not enough to escape the halo, CR pressure gradients can drastically change their dynamics as winds escape from the galaxy.
This section has described how CRs affect outflows, providing further insight on the consistently higher outflow rates that the CRMHD run exhibits when compared to the HD and MHD runs.Here we perform a more detailed analysis of a particular feedback event in the CRMHD simulation to understand how CRs are locally driving the observed differences.In Fig. 11 we show edge-on density-weighted projections of nut during a feedback event at  = 2.41 ( ∼ 2.9 Gyr).This time was selected because of the presence of a corresponding peak in the outflow rate (i.e.∼ 30 M ⊙ yr −1 , see Fig. 7) and the CRsupported 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  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 g cm −3 in black (red) and at 10 −22 g cm −3 in red (blue).These panels uphold the canonical understanding of wind launching in disk galaxies: at the centre lies a cold and dense disk from which high velocity, hot outflows carve biconical perpendicular channels through the warm ionised CGM.This particular outflow is preceded by previous plumes of hot gas visible at larger radii.There are also regions of outflowing gas that do not have a significantly higher temperature than the surrounding CGM (see outflowing region in the top left corner of the  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. −∇    =   , for thermal ( ther ) and CR ( CR ) pressures.To help the analysis, we keep the density contours (green and blue lines) of the temperature projection in these panels.We first examine the hot outflows (top left panels), which display in both thermal and CR pressure support maps a clear biconical morphology with opening angle of ∼ 30 • .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 disk at different times.As such, the front of a previous bubble appears as a positive gradient in the radial direction.These features are also present in the pressure support maps of the HD and MHD simulations (not shown).This reaffirms the selection of the positive support part of these outflows (see Fig. 10) as a method to identify galactic winds.In agreement with the peak observed at  ∼ 2.9 Gyr in Fig. 10, the thermal and CR supports are much larger than unity, both reaching similar magnitudes.Nevertheless, large positive or negative thermal and CR pressure gradients do not always spatially coincide.We surmise that such offsets could be due to CR streaming, as on ∼ kpc scales it can be an important driver of CR transport and displace CRs with respect to the gas, while maintaining the CR pressure gradient (e.g.Dubois et al. 2019).
The warm ionised phase of the gas is the dominant volume-filling phase of the outflow, and has a more spherically isotropic morphology.This phase does not display strong alternating features from consecutive wind bubbles, but rather a smooth pressure gradient.This pressure support extends to large radii and is dominated by CR pressure.In fact, the warm ionised, warm neutral and cold phases exhibit the same support morphology, which suggests that they are all sourced by a common underlying CR pressurised halo that surrounds the galaxy.However, based on the mild support of CR pressure on the warm ionised phase in Fig. 10, we claim that this phase will probably rain down as a galactic fountain given the low average radial velocity (see top right panel of Fig. 9).The four bottom panels examine the warm neutral (left) and cold (right) outflows.Gas in these phases is distributed in small structures that resemble cooling filaments.Blue contours are included to indicate the limits of the hot outflow cones.These colder phases appear to also be moving perpendicularly to the disk, experiencing a uniform CR pressure gradient that dominates their gravitational acceleration.This indicates that these cold regions of gas are still efficiently accelerated by the surrounding hotter phase, and hence entrained in the wind.The radial distribution of these outflowing phases reveals that they decrease in fraction of volume occupied as a function of radial distance, thus allowing to better showcase their clumpy nature.These small outflowing gas structures do traverse the 0.2  vir,DM shell we employ to measure outflows (illustrated by a dashed purple circle) in a timescale shorter than our snapshot time resolution, naturally resulting in the fluctuating behaviour of warm neutral and cold phases present in Figs. 9 and 10.

CONCLUSIONS
In this work, we have performed a detailed study of the effects of CR feedback on a Milky Way-like galaxy using a suite of new zoomin cosmological simulations.The main galaxy in our zoom region, nut, was simulated down to  = 1.5 with the ramses code.Using adaptive mesh refinement, we simulated the multi-phase ISM with a spatial resolution of ∼ 23 pc per full cell size.Our hydrodynamical simulation (HD) includes a model for a resolved ISM, with a magneto-thermo-turbulent (MTT) star formation prescription and a 'mechanical' SN feedback of a canonical strength.In our magnetohydrodynamic simulation (MHD) we additionally followed the magnetic fields with the constrained transport method which are amplified from a primordial seed of 3 × 10 −12 G. Our cosmic ray MHD simulation (CRMHD) further models CRs injected by supernova (SN) events with an efficiency of 10%.In this first work in a series of papers analysing the role of CRs in the nut suite, we focused on the impact of our full CR physics which includes advection, losses, anisotropic diffusion and streaming on the ISM and the launching of outflows.Our main results can be summarised as follows.
(i) Including CRs significantly reduces the resulting stellar mass of nut compared to the no-CR models, yielding good agreement with abundance matching models (e.g.Behroozi et al. 2013;Moster et al. 2013) without the need for a calibrated boost to the SN energy deposition or calibrated galactic winds.This reduction in stellar mass is significant and it reaches approximately an order of magnitude at high redshift  ≳ 3 before tapering to a factor ∼ 4 at  ∼ 1.5 (Section 3.1 and Fig. 4).
(ii) CRs pressure gradients support cold, low density gas against collapse, making the warm phase significantly more dominant in the ISM (Section 3.2).CRs also deplete an important fraction of the mass locked in the cold, dense phase (i.e. ≳ 10 −21 g cm −3 and  < 10 3 K) present in non-CRs simulations, and hence significantly decrease the amount of gas efficiently forming stars.In fact, galactic gas is rarely found below 10 −27 g cm −3 and the ISM has a smoother gas distribution within the disk.
(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 environments where SN events take place are modified.In the CRMHD simulation, SN explosion sites have a bi-modal distribution that both peaks at and reaches considerably lower densities than in the HD and MHD runs (see Fig. 6).This reduces the amount of SN energy lost due to cooling, increases the efficiency of this feedback mechanism and drives stronger outflows.
(iv) The mass outflow rates in the CRMHD simulation are significantly larger at all cosmic times with respect to the no-CR runs, specifically with a higher mass loading factor  especially at high redshift.In contrast to the non-CR simulations where the dominant phase of outflowing gas mass is hot (i.e. > 10 5 K), CRMHD outflows are dominated in mass by the warm ionised phase (i.e.9000 <  < 10 5 K) and have a significant contribution from the cold phase (see Section 3.3.2and 3.3.3).Furthermore, hot outflows are biconical in structure whereas the colder phases embedded in the hot outflow are filamentary and clumpy.This is a direct prediction for upcoming JWST observations of multi-phase galactic outflows expected soon.
(v) Outflows in the CRMHD simulation experience further acceleration due to the CR pressure gradient, with a secondary contribution from the thermal pressure.Outflows in the warm and cold phases have a significant support from the CRs, while the thermal pressure support remains negligible.Therefore, outflows experience further acceleration within the CGM of our simulated galaxy (see Section 3.3.4and Fig. 11).
We find CRs to have a large impact across various properties of the simulated galaxy and its ISM.While we include multiple important CR physical components in our simulations, the specific impact of each is not fully explored.A review of these roles is left to the next paper in this series, which will analyse the effects of this CR feedback on the thermodynamic state of the CGM outside the nut galaxy (Rodríguez Montero et al. in prep.).
While simulations, such as the ones presented here, are revealing CRs to be a central agent in galaxy formation, our CR models are far from being complete and our results suffer from multiple caveats, both physical and numerical.For instance, we still lack a self-consistent coupling of CR transport with the local thermodynamic state of the gas (e.g.Armillotta et al. 2022).Furthermore, CR physics may well be affected by more realistic and complete gas cooling prescriptions (Nuñez-Castiñeyra et al. 2022;Katz et al. 2022) and different refinement criteria in the CGM (Bennett & Sĳacki 2020;Rey et al. 2023).Milky Way-like halos at high redshift can have their star formation strongly altered by other forms of feedback such as that from stellar winds (Agertz et al. 2021), AGN (Koudmani et al. 2022), or radiation from young stars (Rosdahl & Teyssier 2015;Hopkins et al. 2020;Emerick et al. 2018).Their self-consistent interaction with CR feedback remains a topic for further investigation (e.g.Farcy et al. 2022;Martin-Alvarez et al. 2022).
Despite these caveats, our work highlights the fundamental role CRs likely play in the formation and evolution of Milky Way-like galaxies.In particular, they seem to be a central contributor to the feedback loop regulating star formation and driving galactic outflows, providing a promising avenue to form galaxies with more realistic physical properties without having to resort to ad-hoc feedback model calibrations.

Figure 2 .
Figure2.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  = 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  50 obtained by fitting the RGB images with a 2D Sérsic profile.For each simulation the stellar mass within 0.2  vir,DM and galaxy compactness log Σ 1.5(Barro et al. 2013), is listed, with a larger value indicating a more compact stellar emission.Compared to HD and MHD simulations, the disk of nut shows much more structure in the CRMHD simulation, dominated by a more diffuse ISM with less compact and dense clumps.The CGM is also strongly altered, with the cold dense filaments of HD and MHD replaced by a diffuse and cold halo more uniformly enriched by metals.The mock observations of HD and MHD show a galaxy dominated by an old stellar population and a large bulge, while the nut galaxy in CRMHD is transformed into a less compact system dominated by a disk.In this case, star forming regions are scattered across spiral arms and clumps all over the disk.

Figure 3 .
Figure3.Evolution of the stellar (left panel) and gas (right panel) masses contained within the galactic region (i.e.< 0.2 vir,DM ) of the nut galaxy in the HD (blue line), MHD (purple line) and CRMHD (green line) simulations.The right panel distinguishes the total gas mass (solid line) from the cold gas mass (dashed line), as defined by the constant entropy limit of the cold phase of the ISM byGent (2012).In all runs stellar mass grows rapidly up to  ∼ 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 ∼ 1 dex at  ≳ 4. Below  ∼ 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  = 1.5 to be similar to that in the HD and MHD runs, whilst the stellar mass remains 4 times lower.

Figure 4 .
Figure 4. Stellar mass to halo mass function of the nut galaxy for our fiducial models: HD (blue), MHD (magenta) and CRMHD (green) shown in 4 redshift ranges ( ∼ 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) (blueshaded region) and Behroozi et al. (2013) (orange-shaded region), evaluated at  = 1.5.Coloured dashed orange and blue lines indicate the stellar masses at which observational mass functions are subject to incompleteness.The HD and MHD runs are able transform the majority of their baryonic content into stars, while the CRMHD run from early on ( ∼ 8) sits ∼ 1 dex below, approaching the prediction by Behroozi et al. (2013) at  ∼ 1.5.

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

Figure 6 .
Figure6.Distribution of gas densities in which SN events take place in the redshift range  = 6.0 − 3.0 (top panel) and  = 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 ∼ 10 −27 g cm −3 and a peak at ∼ 10 −23 g cm −3 .SN events taking place at lower densities during this period are more effective in suppressing star formation.This bimodality for CRMHD is still present in the lower redshift range (bottom panel), but the contribution of the low density peak is significantly suppressed, with the overall distribution resembling more the one obtained for the HD and MHD simulations.

Figure 7 .
Figure 7. Top panel: Outflow rate measured in a thin spherical shell at 0.2  vir,DM for the HD, MHD and CRMHD runs.The histograms show the instantaneous mass-flow rate through the thin shell, while the solid lines are the running median of these measurements using a filter of size 200 Myr.Middle panel: Star formation rate (SFR) averaged over time bins of 100 Myr.Bottom panel: Mass-loading factor  computed using the running median results of  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 all cosmic time explored.Given the large suppression of star formation for  ≥ 3, the CRMHD simulation has outflow mass-loading consistently exceeding unity, only reached by the HD and MHD simulations during the merger-induced starburst at  ∼ 3 − 2.

Figure 8 .
Figure 8. Histogram of mass outflow rates as a function of gas radial velocity measured in a thin shell at 0.2  vir,DM and ⟨⟩ = 3.94 (top panel) and ⟨⟩ = 2.46 (bottom panel).We compare the HD (blue), MHD (purple),

Figure 9 .
Figure9.Outflow budget and kinematics of hot (top left pair), warm ionised (top right pair), warm neutral (bottom left pair) and cold (bottom right pair) phases for the HD (blue), MHD (purple) and CRMHD (green) runs measured at 0.2  vir,DM .Each pair of panels show the instantaneous outflow rate (top) and the  outflow -weighted radial velocity of outflowing gas (bottom).The  outflow panels show the same quantities as Fig.7: low opacity histogram is the instantaneous value, while the solid lines show the running median with a filter of 200 Myr.In the   panels, the escape velocity of each model at the shell radius is showed with dashed lines, for ease of comparison to the median  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 −1 for  ≥ 3 and ∼ 70 km s −1 for  ≲ 3, dominating the outflow rate in HD and MHD simulations.The hot phase in CRMHD reaches the largest outflow velocities (∼ 150 − 200 km s −1 for  ≥ 2.5 and ∼ 100 km s −1 for  ≲ 2.5), with a non-negligible fraction of the outflowing gas mass exceeding  esc consistently up to  ≥ 2.5.The dominant mass-flow rate is coming from the warm ionised phase in the CRMHD, while above  ∼ 5, warm neutral and cold phases largely contribute to the total outflow rate.Warm neutral and cold phases are scarce in the HD and MHD simulations, usually identified as low velocity peaks which can matched with crossing gas features of orbiting satellites.However, in the CRMHD run these outflows while slower than the hot gas can still reach ∼ 100 km s −1 at lower redshifts.
Figure10.Pressure gradient support against gravity for outflowing gas separated into different temperature phases for HD (blue), MHD (purple) and CRMHD (green) simulations.The dividing line between the shaded and unshaded region indicates marginal support.In the shaded region, gravity dominates the dynamics, while in the unshaded region, gas is accelerated outwards by the pressure gradients.Outflowing gas is measured in a thin shell at 0.2  vir,DM and divided into the phases described in Section 3.3.2according to temperature.Shaded, coloured bands represent the second and fourth quartiles of the underlying distribution.The left column corresponds to the hot gas, with each row corresponding to the HD (top), MHD (middle), and CRMHD (bottom) simulations, respectively.While below  ∼ 3 the thermal pressure struggles to accelerate gas in the HD and MHD simulations, CR pressure consistently supports and accelerates the hot gas across all cosmic time.Only CR pressure can mildly support the warm ionised phase between 3 <  < 6, but it supports and accelerates warm neutral and cold gas for  ≤ 5.