Spectrally resolved cosmic rays - III. Dynamical impact and properties of the circumgalactic medium

Cosmic rays (CRs) are dynamically important in the evolution of galaxies by regulating star formation and powering galactic outflows. However, to what extent CRs regulate galaxy formation depends on the coupling strength of CRs with the ambient plasma and the effective CR transport speed. Moreover, both properties sensitively depend on the CR momentum, which is largely unexplored in three-dimensional hydrodynamical simulations. We perform magneto-hydrodynamical simulations of entire galaxies with masses ranging from 10 10 to 10 12 M ⊙ and compare dynamically coupled CRs in the grey approximation with a spectrally resolved model that includes CR momenta from 0 . 1 GeV c − 1 to 100 TeV c − 1 . We find that hadronic cooling of CRs dominates over Alfv´en cooling, with the latter emulating CR losses as a result of streaming of CRs down their pressure gradient. While star formation rates and galaxy morphologies are only mildly affected by the spectral CR modelling, mass loading factors of galactic outflows can differ by up to a factor of four in dwarf galaxies. All simulated low-mass halos ( M = 10 10 , 10 11 , and 3 × 10 11 M ⊙ ) drive strong outflows, where CR transport is temporally dominated by advection. In contrast, the Milky Way-mass galaxy with M = 10 12 M ⊙ does not drive sustained outflows, so that CR transport is entirely dominated by diffusion. The effective energy weighted diffusion coefficients vary by two orders of magnitude from the canonical energy-weighted values of ⟨ D ⟩ e cr ∼ 10 28 cm 2 s − 1 in the disc up to 3 × 10 29 cm 2 s − 1 in the circumgalactic medium, where we observe substantial temperature and CR pressure differences between our grey and spectral CR models.


INTRODUCTION
CRs constitute an important energy reservoir in the interstellar medium of star forming galaxies (Strong et al. 2007;Grenier et al. 2015;Zweibel 2017).Their energy density is comparable to that in the other components (i.e., the thermal, magnetic, and turbulent energy densities) of the interstellar medium (ISM, Cox 2005;Naab & Ostriker 2017), which enables them to power dynamically relevant processes (see Ruszkowski & Pfrommer 2023, for a review).CRs are known to regulate star formation in the densest regions of the ISM via CR heating and ionization (e.g., Ivlev et al. 2018;Phan et al. 2018;Padovani et al. 2020) and drive powerful winds from the ISM into the halo (e.g., Naab & Ostriker 2017;Veilleux et al. 2020).They also constitute a significant fraction of the total pressure in the circum-galactic medium (CGM, Salem et al. 2016;Buck et al. 2020;Ji et al. 2020).
CRs have also been included in smaller scale simulations of the ISM, which explicitly resolve the multi-phase nature of the gas and resolve the dynamical launching of the outflow from the disc more accurately than the galactic models (Girichidis et al. 2016;Simpson et al. 2016;Girichidis et al. 2018;Farber et al. 2018;Commerçon et al. 2019;Armillotta et al. 2021;Rathjen et al. 2021;Simpson et al. 2022).These small-scale models reinforce the conclusion that CRs can have an important dynamical impact on the ISM and are crucial for launching and accelerating galactic outflows.
Whether CRs are most effective in driving the warm or hot phase of the ISM, and for which galaxies at which cosmic time CRs have their largest impact, sensitively depends on the CR transport mechanisms and their coupling with the ambient plasma, namely gyro-resonant CR-wave interactions.While CRs below ∼ 100 GeV are expected to scatter off resonant MHD waves generated predominantly by the gyroresonant (streaming) instability (Kulsrud & Pearce 1969) or the intermediate-scale instability (Shalaby et al. 2021(Shalaby et al. , 2022(Shalaby et al. , 2023)), CRs of higher energies are expected to scatter off MHD waves driven primarily by extrinsic turbulence (Aloisio et al. 2015;Zweibel 2017;Evoli et al. 2018).In addition to the gyro-resonant (streaming) instability, CRs can amplify magnetic fields via the non-resonant instability (Bell 2004), which can then reduce the effective transport of CRs on scales close to their sites of acceleration (which are unresolved in the present work).Describing CR transport by a two-moment formulation of CR hydrodynamics enables a self-consistent computation of the macroscopic CR transport speed and the associated diffusion coefficient, using either a steady-state approximation for the resonant wave energy (e.g., Jiang & Oh 2018;Hopkins et al. 2021), or a formulation that dynamically follows the resonant Alfvén wave energy in response to the gyroresonant CR streaming instability and various collisionless damping processes (Thomas & Pfrommer 2019, 2022;Thomas et al. 2020Thomas et al. , 2021)).
Furthermore, most CR cooling and transport effects are strong functions of the CR energy.Because the distribution of CRs follows a power law energy spectrum ranging from MeV to 10 20 eV, their interactions with the ambient plasma vary strongly depending on the energy regime.Low-energy CRs (≲ 0.1 GeV) have an increased cross section with the thermal gas and provide most of the heating in dense star-forming cores (e.g., Girichidis et al. 2020;Padovani et al. 2020).CRs in the GeV energy range cool less efficiently.As a result, they are primarily responsible for carrying the pressure of the entire population for realistic values of the momentum spectral index of 2.2 to 2.8.Consequently, they are responsible for driving galactic outflows.While the high-energy part of the CR spectrum with energies ≳ TeV is typically subdominant in total energy and therefore dynamically not relevant, these high-energy CRs produce γ rays, which is of prime interest for imaging air and water Cerenkov telescopes such as CTA, HAWK, and LHASSO.
In particular, CR ions interact hadronically with the ambient medium and produce pions that decay into γ rays, neutrinos, and secondary electrons/positrons.Similarly to primary, shock accelerated electrons, these secondary leptons produce radio synchrotron emission and inverse Compton γ-ray emission, which enables studying the underlying CR populations via their non-thermal multi-frequency spectra from the radio to PeV γ-ray range (Kotera & Olinto 2011;Werhahn et al. 2021a).A testimony of this is the far-infrared-γ-ray relation at GeV energies (Ackermann et al. 2012;Rojas-Bravo & Araya 2016;Ajello et al. 2020) that has been modelled in three-dimensional MHD simulations (Pfrommer et al. 2017b;Chan et al. 2019;Werhahn et al. 2021b;Buck et al. 2020;Nuñez-Castiñeyra et al. 2022), indicating that it allows probing the CR ion population via the dominant pion-decay γ-ray emission, and indirectly its ability to provide dynamical feedback through the calorimetric fraction (i.e., the energy fraction that is radiated away via hadronic CR interactions and thus not available any more for feedback).Similarly interesting is the far-infrared-radio relation (van der Kruit 1971;Bell 2003), which enables probing CR electron calorimetry as well as transport and ageing processes (Thompson et al. 2006;Lacki et al. 2010;Werhahn et al. 2021c;Pfrommer et al. 2022).
Most previous hydrodynamical models only include CRs in the grey approximation, i.e., they follow the momentumintegrated total CR energy.The transport coefficients as well as the cooling processes are treated in an effective manner without any momentum dependence while assuming a globally universal steady-state CR spectrum (e.g., Pfrommer et al. 2017a).However, numerical simulations that include a full CR spectrum (Girichidis et al. 2020, hereafter Paper I) suggest that the shape of the spectra differ perceptibly in different regions of the galaxy as well as over time (Girichidis et al. 2022, hereafter Paper II).This finding is also supported by a spectral analysis in post-processing of CR-MHD models (Werhahn et al. 2021a).
In this study, we extend the previous models of spectrally resolved, dynamically coupled CRs (Paper I; Paper II) in galaxy simulations using the moving mesh code Arepo (Springel 2010;Weinberger et al. 2020).We compare spectrally resolved CRs to their grey counterparts and investigate the dynamical and thermal impact of CRs in different halos with masses ranging from dwarfs (M = 10 10 M⊙) to Milky Way analogues (M = 10 12 M⊙).We emphasize that the spectral model in (Paper II) erroneously did not account for CR cooling.The old model therefore marks an extreme case of energy dependent CR transport, which we extend here to also include CR cooling.
The outline of the paper is as follows.In Section 2 we briefly summarize the basics of the spectral model (Paper I; Paper II). Section 3 describes the galaxy setups and details the simulations performed.An analysis of the cooling processes is presented in Section 4, followed by an analysis of the star formation rate and the outflows in Section 5. We discuss the effective diffusion coefficients and the dominant transport mechanisms in Sections 6 and 7, respectively.The time averaged accelerations and their connection to outflows are discussed in Section 8. Finally, we explore CGM properties in Section 9, before finishing the paper with a discussion (Section 10) and conclusions (Section 11).In Appendix A, we show the effective adiabatic index and the associated pressure of the CR distribution and estimate the error on the gas dynamics associated with our simplified Alfvén cooling approach in Appendix B. In Appendix C, we show the morphology of our dynamo-grown magnetic field.

Spectral CR description and coupling to MHD
We use the spectrally resolved CR-MHD model developed in Paper I (spectral scheme) and Paper II (coupling to MHD), of which we highlight the main concept and features here.
Table 1.List of simulations.Shown are the simulation name, the halo mass M , and whether the simulation includes spectral CRs (in this case, N bins denotes the number of spectral bins).The next two columns list the diffusion coefficient along the magnetic field, which corresponds to D(1 GeV c −1 ) in the spectral case, and the momentum scaling index that is defined via D ∝ p δ .The last column indicates whether Alfvén cooling is included or not.
no. name mass spectral M 10 10 -grey4×10 other losses and Fermi II acceleration where f = f (x , p, t) = d 6 N/(dx 3 dp 3 ) denotes the isotropic part of the CR distribution function, with the spatial coordinates x , the momentum p = |p|, and the time t.The mean velocity of the gas is v, the spatial diffusion tensor is Dxx = Dxx(x , p, t), and the diffusion coefficient in momentum space is Dpp = Dpp(x , p, t).CR losses and sources are expressed as b l = b l (x , p, t) = dp/dt and j = j(x , p, t), respectively.The losses b l include Coulomb and hadronic in-teractions.Since we neglect diffusion in momentum space (Dpp = 0), we will omit the subscript of the spatial diffusion coefficient in the following.The momentum space is discretised using Nspec logarithmic bins with bin face momenta p i−1/2 and bin centred counterparts pi = √ p i−1/2 p i+1/2 .The particle distribution function f is described by piecewise power laws where p i−1/2 is the momentum and f i−1/2 is the corresponding amplitude at the left-hand face of momentum bin i.The slope in bin i is denoted by qi and θ is the Heaviside function.We note that the spectrum does not need to be continuous at the bin faces, which increases numerical accuracy and avoids instabilities (see Paper I). Solving for the two degrees of freedom per bin (f i−1/2 and qi) requires solving for two independent quantities related to f .We choose the number and energy density, where is the kinetic energy of individual protons and c is the speed of light.The CRs are coupled to the MHD equations via the gas velocity v for advection, and the total CR pressure for the dynamical coupling in the Riemann solver (see Pfrommer et al. 2017a).

Spatial diffusion, streaming approximation and cooling
We use momentum dependent anisotropic diffusion along the direction of the magnetic field using the implicit solver described in Pakmor et al. (2016b).An improvement in terms of spectral stability (Paper II) allows us to only diffuse the CR energy density -rather than both number and energy density -and reconstruct the spectrum in a stable and energy conserving manner after the diffusion step.Diffusion is applied to each spectral bin with a bin-centred diffusion coefficient parallel to the magnetic field1 Diffusion itself is an energy conserving process, and in order to emulate CR streaming along the magnetic field, we additionally account for CR Alfvén wave losses.The streaming instability effectively reduces the mean CR transport speed, vst, to the Alfvén speed, where B denotes the magnetic field and ρ is the gas mass density (assuming full ionisation We apply ΛA to the total CR energy after the diffusion step and rescale the spectrum to match the new CR energy.In Appendix B, we show our simulated CR spectra and demonstrate that the error on the gas dynamics associated with this rescaling is expected to be small.In addition, we account for Coulomb, ionization and hadronic losses as described in Paper I.

GALAXY SIMULATION SETUP
We use the Arepo code, which employs an unstructured moving mesh to decompose the computational domain.The hydrodynamical solver is second-order accurate in space and time (Springel 2010;Pakmor et al. 2016a;Weinberger et al. 2020).CRs are dynamically coupled to the evolution of the gas as described in Pfrommer et al. (2017a) using the advection-diffusion approximation.The spectral solver for the CRs is coupled to the MHD and the diffusion solver as outlined in Paper II. Magnetic fields are included in the ideal MHD approximation (Pakmor et al. 2011;Pakmor & Springel 2013) using the Powell 8-wave scheme for divergence cleaning of the magnetic field (Powell et al. 1999).This enables us to resolve magnetic field growth via a small-scale dynamo already for turbulent injection length scales of order 1 kpc (Pakmor et al. 2017;Pfrommer et al. 2022).CR diffusion uses the semi-implicit diffusion solver by Pakmor et al. (2016c).
For each galaxy we simulate, the numerical setup consists of an isolated collapsing gas cloud that is embedded in a dark matter halo, similar to the setups in Pakmor et al. (2016c) and Pfrommer et al. (2017a).We adopt a baryon mass fraction of Ω b /Ωm = 0.155 (Planck Collaboration et al. 2020) for all our models and vary the total mass of the system.Our four different haloes have total masses of 10 10 M⊙, 10 11 M⊙, 3 × 10 11 M⊙ and 10 12 M⊙ (our Milky Way-like halo) with corresponding initial gas masses of 1.55 × 10 9 M⊙, 1.55 × 10 10 M⊙, 4.65 × 10 10 M⊙ and 1.55 × 10 11 M⊙, respectively.The dark matter and the gas initially follow an NFW profile (Navarro et al. 1997) with a concentration parameter of c200 = r200/rs = 7.Here, rs denotes the characteristic scale radius of the NFW profile.The radius r200 encloses a mass with an average density 200 times the critical density of the universe.While we adopt an external fixed dark matter potential, the gas follows initially a solid body rotation with a dimensionless spin parameter λ = J|E| 1/2 /(GM Here, J is the angular momentum, |E| the total halo energy and G Newton's constant.The mass inside r200 is given by M200, which we abbreviate with M in our setups to keep the simulation names short.We set an initial seed field along the x-axis with a field strength of 10 −4 µG.The halo is sampled with 10 6 gas cells with a uniform target mass resolution of 10 −6 times the total gas mass, which ranges from 1.55 × 103 M⊙ for the galaxy with the lowest mass up to 1.55 × 10 5 M⊙ for our Milky Way-like system, respectively.Adjacent cells are only allowed to differ by a factor of ten in volume. We include radiative cooling and star formation employing the pressurized ISM model by Springel & Hernquist (2003), which includes a stochastic star formation prescription that reproduces the observational relation between the gas surface density and the star formation rate surface density (Kennicutt 1998).The subgrid ISM model assumes the hot and cold phase of the ISM to be in pressure equilibrium.The total ISM pressure is modeled with a stiff equation of state above a threshold gas density for star formation, nSF = 0.13 cm −3 .Star formation is implemented as an injection of a star particle with instantaneous, implicit stellar feedback as the effective equation of state used as subgrid ISM model does not require modelling explicit thermal/kinetic feedback.We inject CRs with an efficiency of 10 per cent of the canonical supernova (SN) injection energy of 10 51 erg (see, e.g., Helder et al. 2012;Morlino & Caprioli 2012;Ackermann et al. 2013) with an injection spectrum of f (p) ∝ p 4.2 .However, detailed comparisons of numerical models with observations of individual SNe suggest that also lower efficiencies of only 5 per cent are plausible (Pais et al. 2018).Initially, there are no CRs in the models, so all CR energy is injected by SNe.
We present results from 32 simulations, which are listed in Table 1.We distinguish between grey CR models similar to the galaxies in Pakmor et al. (2016c); Pfrommer et al. (2017a) using different diffusion coefficients ("grey"), spectral CR models ("spec") and control runs without CRs ("noCR").We further differentiate between models with and without Alfvén cooling (ΛA and ΛA = 0).For the Milky Way-like galaxies, we furthermore study two different momentum scalings of the diffusion coefficient (Eq.( 7)) with δ = 0.3 to 0.5 (see Evoli et al. 2020;Werhahn et al. 2021a,b,c, for a discussion of the momentum scaling).The momenta for the spectral bin edges, the bin centres, and the corresponding bin-centred diffusion coefficients are listed in Table 2.
For the spectral models we use 12 momentum bins ranging from 0.1 GeV c −1 to 100 TeV c −1 .We inject CRs with a power-law source spectrum, f (p) ∝ p −4.2 , into the local environment of newly created star particles and use a spherical top hat filter containing the closest 32 mesh cells surrounding the star.We run all simulations for a total time of 3 Gyr.

CR COOLING PROCESSES
Before analyzing the simulations we can estimate the expected scalings of the CR energy distribution and the cooling times with gas density.For adiabatic changes of the CRs, we expect the CR energy density and pressure to scale as ecr, Pcr ∝ ρ 4/3 , where we adopt the adiabatic index γ = 4/3 for a relativistic fluid. 2 At high gas densities, non-adiabatic CR interactions with the gas as well as active CR transport processes such as diffusion are expected to cause a deviation of this adiabatic scaling (see, e.g., figure 16 of Girichidis et al. 2018).For hadronic interactions with ambient nuclei of the ISM, the loss rate of kinetic CR energy reads where c denotes the speed of light, n the target nucleon density of the ISM, σpp is the total pion cross section and Kp ≈ 1/2 represents the inelasticity of the reaction (Mannheim & Schlickeiser 1994).The threshold momentum above which pions can be produced hadronically is p thr ≈ 0.78 GeV c −1 (θ is the Heaviside step function).The cooling rate thus scales linearly with n, and the cooling time T / Ṫ is independent of the CR momentum 3 with Thus, the hadronic cooling rate has the following scaling with gas density for adiabatic CRs: In order to derive the scaling of the Alfvén cooling rate with gas density, we investigate the individual terms that enter equation ( 9).Assuming magnetic flux-freezing for isotropic changes of the volume, we have |B| ∝ ρ 2/3 in the kinematic limit (e.g., Spruit 2013), so that |vA| ∝ ρ 1/6 .We approximate the gradient of the CR pressure as ∇Pcr ∼ Pcr/Lcr, where Lcr is the gradient length of the CR pressure, and assume Lcr ∝ ρ −1/3 for the scaling of the CR gradient length.Assuming adiabatic changes of the CR pressure, the Alfvén wave losses are expected to scale as and the corresponding cooling time scales for adiabatic CRs as which shows a weaker density scaling in comparison to the hadronic CR interactions (cf.Eq. 11).
Figure 1 illustrates from left to right the distribution of the CR energy density, the hadronic as well as the Alfvén cooling rate as a function of density for simulation "M 10 11 -spec-ΛA" averaged from t = 1 − 2 Gyr.Averaging over different time intervals does not yield any differences (except for the very beginning of the simulation when the magnetic dynamo has not yet saturated, Pfrommer et al. 2022).Furthermore, the other simulations behave similarly.Colour coded is the fraction of the total CR energy. 4The dashed lines show the median along the ordinate, the solid lines indicate the analytic Figure 2. Illustration of the cooling times for simulation "M 10 11 -spec-Λ A " at t = 1.5 Gyr.From left to right we show gas density, the total CR energy density, and the cooling times for Alfvén cooling and hadronic cooling, respectively.The Alfvén cooling time is significantly larger than the hadronic counterpart in most of the volume.Only in the centre, t cool,A is comparable to the dynamical time.This is contrary to the hadronic case, where most of the disc has cooling times perceptibly shorter than the dynamical time, which is a necessary condition for CR calorimetry.
scaling of the corresponding variable with ρ while assuming adiabatic CRs.The scaling of the energy density closely follows the adiabatic approximation for low densities.Above ρ ∼ 10 −26 g cm −3 the scaling starts to flatten, which is related to non-adiabatic effects such as non-adiabatic cooling in the dense gas and diffusion of CRs out of the dense gas, where they have been injected.The relatively tight distribution of ecr around the median leads to a narrow distribution each snapshot has an equal weight in the averaged distribution irrespective of its total CR energy content.We then compute the 2D histogram and divide it by the number of snapshots.
of the hadronic cooling rate, which follows the adiabatic CR scaling up to that density, at which point the adiabatic approximation for ecr breaks down.
The comparison of the hadronic with the Alfvén cooling rate in the right-hand panel illustrates two important features.First, the distribution shows a significantly larger scatter, which is the result of larger local differences in the magnetic field strength and local CR pressure gradients.Second, it shows that for densities above ρ ∼ 10 −28 g cm −3 , the hadronic cooling dominates over the Alfvén cooling rate.In the central region of the galaxy at densities above ρ ≳ 10 −25 g cm −3 the median of the hadronic cooling rate is ap-  legend).The differences between the models at the identical halo mass is negligible.The different halo masses mainly occupy different ranges in the gas column density.Below a surface density of Σgas ≈ 10 M ⊙ pc −2 the star formation rate reduces compared to the classical relation, which is a well-known observed feature (e.g., Bigiel et al. 2008;Leroy et al. 2008) proximately two orders of magnitude larger than the Alfvénic counterpart, which explains the small differences between the simulations with and without Alfvén cooling.
We illustrate the spatial distribution of the Alfvén and hadronic cooling times in Fig. 2, where we show the gas density, the CR energy density, the Alfvén cooling time and the hadronic cooling time in cuts through the centre of the disc for the simulation "M 10 11 -spec-ΛA" at t = 1.5 Gyr.Both the density and the CR energy density show a relatively smooth distribution across three orders of magnitude.The measured cooling times t cool,A range from ∼ 100 Myr in the centre to values larger than the Hubble time in the outskirts of the galactic disc.Only in the innermost ∼ 5 kpc the cooling times are comparable or shorter than the dynamcial time of the disc (≲ 100 Myr).For the hadronic counterpart (t cool,hadr ) we find cooling times ranging from ∼ 10 Myr to ∼ 10 Gyr.For most of the volume, t cool,hadr ≪ t cool,A .Furthermore, for most parts of the disc, the hadronic cooling time is shorter or comparable to the dynamical time of the disc, which is a necessary condition for CR calorimetry.

STAR FORMATION AND OUTFLOW RATE
The dynamical evolution proceeds very similarly for all setups5 .The initial gas cloud rapidly starts cooling inside out, collapses and forms stars in close proximity to the centre of the halo.The strength of this starburst phase depends on the mass of the system, with higher star formation rates for larger masses.The initial rotation of the gas cloud and the conservation of specific angular momentum during the collapse causes the formation of a disc inside out.This first phase of the gravitational collapse is dominated by adiabatic compression of the gas (Pfrommer et al. 2022).The impact of the CRs is limited during this phase since they are only injected during the initial burst of star formation and need time to build up relevant energy densities.
Before investigating the details we confirm that the star formation rate in the simulations scales with gas surface density as expected (Fig. 3).Included are all models at each simulated halo mass over the full time evolution.The differences between the individual models are small (see below).Overplotted is the classical Kennicutt-Schmidt relation (Schmidt 1959;Kennicutt 1998).We identify the newly formed stars in patches of (2 kpc) 2 in the xy-plane and compute the corresponding column density along the z-direction in the range of ±5 kpc around the midplane.The modelled systems with different halo masses mainly occupy different regions in the plot.At low surface densities, we note a significant deviation from the classical relation towards lower star formation rate surface densities, which is broadly in agreement with observations (e.g., Bigiel et al. 2008;Leroy et al. 2008) Figure 4 (left-hand panels) shows the time evolution of the star formation rate for all models.From top to bottom we increase the halo mass.For the early starburst phase we only find strong differences in the star formation rate between the grey and spectral CR models for the M = 10 10 M⊙ halo.The low-mass systems are less strongly dominated by gravitational collapse, so that details in the injection rate, the location and the cooling can alter the temporal evolution.Until a time of t ∼ 0.5 Gyr, the spectral model behaves more similar to the "M 10 10 -grey4 × 10 28 -ΛA" run, whereas afterwards it mimics a decline that is comparable in shape to the "M 10 10 -grey10 28 -ΛA" setup, however, with an offset in amplitude so that it produces the largest star formation rate over the simulated time.At late times we also note higher rates for the grey run that includes Alfvén cooling "M 10 10 -grey10 28 -ΛA" compared to "M 10 10 -grey10 28 -ΛA = 0", which stems from the more efficient CR cooling and the resulting weaker CR support against gravitational collapse.For the M = 10 11 M⊙ halos, the differences in the starburst phase are negligible.At a later stage, the star formation rates in the spectral models do not decline as strongly in comparison to the grey counterparts.Towards the end of the simulation the missing Alfvén cooling (dashed lines) again suppresses star formation in comparison to the counterpart runs that include ΛA.For the two largest halo masses, the impact of the gravitationally driven collapse of the halo gas overlies the weaker dynamical differences between the different CR models.Hence, the temporal evolution of the star formation rate does not indicate strong systematic differences.The systematic reduction of the star formation rate due to missing Alfvén cooling as observed in the smaller mass haloes is subdominant in the temporal scatter of the largest halo.At late times the star formation rate in the grey run with D = 10 29 cm 2 s −1 declines noticeably faster than in all other runs.
We quantify the outflows as a function of time in M 10 10 -grey10 28 -Λ A = 0 M 10 10 -grey10 28 -Λ A M 10 10 -grey4×10 28 -Λ A M 10 10 -grey10 29 -Λ A M 10 10 -grey3×10 29 -Λ A M 10 10 -spec-Λ A Figure 4. Left: Time evolution of the star formation rate for all halo masses (top to bottom).The peak of the star formation rate hardly differs between the grey and spectral models, except for the halo with the lowest mass (M = 10 10 M ⊙ ).At late times, the spectral models show higher star formation rates for the two setups with low halo masses (M = 10 10 M ⊙ and M = 10 11 M ⊙ ) but not for the two more massive ones.Right: Outflow rate as a function of time for all halo masses measured at a height of 10 kpc above and below the midplane through a circular area with radius r = 15 kpc.We measure noticeable outflows in all but the highest mass halo (M = 10 12 M ⊙ ), which is dominated by infall.The strongest outflows are driven by the grey models with diffusion coefficients D = 4 × 10 28 and 10 29 cm 2 s −1 with a peak at around t ∼ 0.5 − 1 Gyr.For the most massive halo (M = 10 12 M ⊙ ) the grey model with the highest diffusion coefficient (D = 3 × 10 29 cm 2 s −1 ) can drive an outflow for a short period of time.phase, i.e. for t ≲ 0.5 Gyr, where all simulations show negative outflow rates.In the case of the low halo mass (10 10 M⊙), the stellar feedback reverts the flow to an effective outflow at t ≈ 0.5 Gyr.The strongest outflows for these dwarf galaxies are driven by model "M 10 10 -grey4 × 10 28 -ΛA" whilst the grey model with a smaller diffusion coefficient, "M 10 10 -grey10 28 -ΛA", and the spectral model "M 10 10 -spec-ΛA" are almost indistinguishable for the entire simulation time.Stronger mass outflow rates are triggered in the 10 11 M⊙ models.Again, the grey model "M 10 11 -grey4 × 10 28 -ΛA" drives the strongest outflows.The inclusion of Alfvén cooling does not alter the outflow properties.The net outflow remains positive for most of the simulated time, but declines to irrelevant net values towards the end of the simulation.For a halo mass of 3 × 10 11 M⊙ all models are able to lift gas above the measurement height of 10 kpc during the initial starburst.Here, the grey models with D = 4 × 10 28 , 10 29 , and 3 × 10 29 cm 2 s −1 drive stronger peak outflows.In the case of the spectral model (M 3×10 11 -spec-ΛA) this short period is followed by weak infall for about 500 Myr.Only after t ≳ 1.5 Gyr, the flow can be reverted to a net outflow.For the highest masses we do not find net outflows for most of the simulated time.Only the grey model with the highest diffusion coefficient (M 10 12 -grey3×10 29 -ΛA) creates a net outflow for a short amount of time after the peak in star formation (∼ 1 Gyr) and for a more evolved stage later on (t > 2 Gyr).
An equally important quantity is the mass loading factor η(t) = Ṁout(t)/ ṀSFR(t).To compute η, we evaluate outflow and star formation rates at the same time despite the fact that the effective stellar feedback that leads to the outflows at a given height has taken place earlier.With typical outflow velocities exceeding 50 km s −1 and a measurement height for the outflows at |z| = 10 kpc, this results in a delay of ≲ 200 Myr, which is significantly smaller than the simulation time.To avoid time scale problems we average the massloading factor between t = 0.5 and t = 2.5 Gyr, which is shown in Fig. 5 as a function of halo mass.We apply a small horizontal shift for better readability.The largest values for η are measured for the smallest halo and η declines to negative values when increasing the system mass from 3 × 10 11 to 10 12 M⊙.The variations between the models span a factor of three to four.For the grey models we note a systematic trend.With increasing diffusion coefficient η first increases and then decreases at a fixed halo mass.Most efficient in driving outflows are the models with D = 4 × 10 28 and 10 29 cm 2 s −1 .An exception is the most massive halo, where the grey model with D = 3 × 10 29 cm 2 s −1 has the largest η.The spectral model and the grey counterpart with a small diffusion coefficient (D = 10 28 cm 2 s −1 ) are comparable.Switching off Alfvén cooling noticeably changes the CR energy content in the centre of the galaxy, where most stars form and where the outflows are primarily launched.As a consequence, the two models with ΛA = 0 show enhanced mass loading factors.The models with a halo mass of 3 × 10 11 M⊙ indicate a similar trend, however, with smaller differences while there is no difference in the most massive halo with 10 12 M⊙.
The formation of an accretion shock that travels from the forming disc into the halo and the accompanying onset of an outflow depends noticeably on the effective diffusion coefficient.We illustrate this behaviour in Fig. 6 for a halo mass of 10 12 M⊙.The simulations with a halo mass of 3 × 10 11 M⊙ show a very similar feature; for the two small halo masses the differences between the grey and spectral runs concerning the accretion shock are very minor, which is a result of the advection dominated CR transport, see below.From top to bottom the rows depict edge-on views of the vertical velocity and the gas density, followed by face-on views of the gas density and the total CR energy density.The four left-hand columns depict the grey models with different diffusion coefficients, and the right-hand column shows the spectral model.In the latter case, the high-energy CRs with diffusion coefficients exceeding most grey counterparts result in faster transport of CR energy into the halo.This triggers a faster expansion of the accretion shock and a faster onset of the outflow.The effective diffusion coefficient of the spectral model at this early stage of the hydrodynamical evolution is comparable to the diffusion in the grey model with D = 10 29 cm 2 s −1 (see next section).We therefore note similarities between all quantities.The faster transfer of CR energy into the outskirts of the galaxy is nicely illustrated in the bottom row of Fig. 6.We note that the overall shape of the galactic disc does not depend on the CR model except for the grey model with the highest diffusion coefficient, which forms a smaller disk and noticeably larger perturbations in the density structure.

EFFECTIVE DIFFUSION COEFFICIENTS
For the spectral models, we can measure the effective diffusion coefficient with which the CR energy is transported.We compute the CR energy weighted value where ecr is the total CR energy.In our discretised spectrum this average corresponds to a sum over the spectral bins, i: The energy weighted diffusion coefficients are shown in Fig. 7 at two different times for different halo masses.The upper panels depict an early evolutionary stage at approximately the peak of star formation (t = 0.2 Gyr) while the bottom panels show a later stage at t = 1 Gyr.For each time we show the edge-on view in the upper row and the face-on view in the lower one.In all cases we use the fiducial spectral model with D0 = 10 28 cm 2 s −1 , in which the diffusion coefficient scales as D(p) ∝ p δ with δ = 0.3.We exclude regions with a total CR energy density (ecr) below 10 −5 eV cm −3 (white patches) because they correspond to pristine conditions in our initial conditions.The white circle indicates the radius enclosing 75 per cent of the stellar mass to highlight the primary regions of CR injection.At the peak of the star formation rate, the central parts of the galaxies are dominated by spectra that closely correspond to the injection spectra with an abundant low-energy CR component.Hence, the energy weighted diffusion coefficient is small with a value ∼ 10 28 cm 2 s −1 .The regions between the star forming sites and the accretion shock From left to right we show cuts through the middle of the box for halo masses of 10 10 , 10 11 , 3 × 10 11 , and 10 12 M ⊙ .We exclude regions that are not yet populated by CRs (white areas).The circles indicate the radius of 75 per cent enclosed stellar mass.At early times we find small effective diffusivities (⟨D⟩e cr ∼ 10 28 cm 2 s −1 ) close to the centre, where new CRs are injected.At larger distances from the centre the spectra are more dominated by high energy CRs, which increases ⟨D⟩e cr to values of ∼ 3 × 10 29 cm 2 s −1 .At late times we note a strong difference between the regions that are dominated by advective transport with low-energy CRs dominating the spectrum and regions of diffusive transport with higher diffusivities.shows values of ⟨D⟩e cr ≳ 10 29 cm 2 s −1 .The relative fraction of the volume with low values of ⟨D⟩e cr is much larger in the lowmass systems in comparison to the higher mass halos so that in model "M 10 10 -spec-ΛA" most of the volume is dominated by diffusion coefficients ⟨D⟩e cr ∼ 10 28 cm 2 s −1 .Contrary, in run "M 10 12 -spec-ΛA" the region with ⟨D⟩e cr ∼ 10 28 cm 2 s −1 is concentrated very closely to the centre.
Low values of ⟨D⟩e cr indicate more abundant low-energy CRs compared to their high-energy counterparts.The existence of large volumes with small effective CR diffusion co-efficients beyond the injection sites (Fig. 7) can have several reasons.First, a likely possibility is advection of CRs, which preserves their spectral shape.A second possibility is adiabatic cooling, which shifts the spectrum to lower momenta and thus gives the low-energy part a larger weight.The large volume with overall lower ⟨D⟩e cr in low-mass halos could also be due to a different calorimetric fraction across the different halo masses.If CRs cool less efficiently at the smaller column densities in dwarf galaxies, low-energy CRs survive for a longer time while high-energy CRs diffuse to larger dis- tances.We note that the three possibilities are connected in a complicated way.A different effective transport is likely to result in different cooling efficiencies.Adiabatic compression or expansion shifts the spectrum and thus changes both the cooling efficiencies as well as the effective transport.We will address the details of the spectral CR shape and the calorimetry of individual energies in a follow-up paper and focus here on the effect of the different transport mechanisms.Highenergy CRs can faster escape the star forming regions, which leads to spectra that are more dominated by high-energy CRs (see Girichidis et al. in prep.).In contrast, low-energy CRs need much longer to escape their acceleration sites.they do not cool as fast as the high-energy ones diffuse away, the spectra peak at higher energies and thus increase the effective diffusivity (see also Werhahn et al. 2023).
At later times (t = 1 Gyr, Fig. 7, bottom panels) the values for ⟨D⟩e cr differ more strongly between the models, which directly reflects their dynamical state.The simulation "M 10 10 -spec-ΛA" drives weak outflows.In spatial regions where CRs are mainly advected and where they have developed a steady state, the CR spectra are dominated by low-to medium-energy CRs, which results in effective diffusivities of ⟨D⟩e cr ∼ 10 28 cm 2 s −1 .This is particularly prominent in the strong outflow of model "M 10 11 -spec-ΛA".The larger gravitational attraction in the halo masses 3 × 10 11 and 10 12 M⊙ results in weaker outflows, and therefore a smaller fraction of advected CR spectra.

CR ADVECTION VS. DIFFUSION DOMINATED SYSTEMS
We investigate the relative importance of advective and diffusive CR transport in the outflow.We note that the role of advection is twofold: while CRs are advected with the fluxfrozen magnetic fields and thus the outflowing gas (in addition to diffusive CR transport along the orientation of the local magnetic field), the observed galactic outflows have been driven by the CR pressure gradients.
For this analysis we focus primarily on the early evolution during the starburst phase.In the case of advection, we only consider radially outflowing velocities, vr > 0, where we define The diffusive speed is computed as where we use the energy weighted effective diffusion coefficient ⟨D⟩e cr from Eq. ( 16). Figure 8 illustrates these different velocities in edge-on views for our fiducial spectral CR model, with increasing halo mass from left to right.From top to bottom we show the gas density, the radially outward pointing component of the velocity vr, the diffusive velocity, and the ratio of the latter two, vr/v diff .In the maps of v diff and vr/v diff we again exclude regions which have not been populated by injected CRs, i.e., regions with ecr < 10 −5 eV cm −3 (white areas).The maps show averages along the line of sight in a slice with a thickness of 5 kpc centred around the centre of the galaxy.For the lowest mass galaxy, the radial velocities indicate the relatively smooth dynamical structure in a coherent outflowing bubble.With increasing halo mass the perturbations in the CGM increase.The two intermediate-mass models show patches of outflowing gas.In the Milky Way-like model, the motions are dominated by infall so that we do not detect large regions of a sustained radial outflow.
The diffusive speeds in the lowest mass halo clearly reflect the coherent radial outflow: the edge to the white areas of negligible CR energy (pristine gas of the initial conditions) approximately coincides with the outer edge of the advective bubble in the map of vr.Furthermore, the diffusive speeds indicate that the outflow is dominated by advection: the locations of the edge of the visible bubble in the vr map correspond to very large diffusive speeds.Those large values arise from strong CR gradients, which are established by advecting most of the CR energy up to this edge.If diffusion were to be dominant throughout the outflow, a larger fraction of the CR energy would be able to diffuse ahead of the radially driven advective outflow, which would lower the CR gradients and the resulting diffusive speeds.This transition is clearly visible for the intermediate-mass models.The Milky Way analogue develops volume filling diffusive speeds between 30 and 100 km s −1 without a clear correlation between the advective and diffusive patterns.This overall behaviour is also reflected in the velocity ratios (bottom row in Fig. 8).The patches of regions where CR transport is primarily advective reduce with increasing halo mass.
We also look at this in a more quantitative manner by computing the volume averaged ratio where i indicates the sum over all cells i with their individual volumes Vi, and V = i Vi is the total volume.We investigate a cylinder with radius 25 kpc and vertical extent ±12.5 kpc centred on the galaxy so that V = π25 2 × 25 kpc 3 . .Vertical profiles averaged over a cylinder with radius r = 15 kpc and height ±30 kpc positioned on the centre of the galaxy at t = 0.3 Gyr (left) and t = 0.6 Gyr (right) for simulation "M 10 11 -spec-Λ A ". From top to bottom we show the acceleration of the gas owing to CRs, their energy density, the gas velocity as well as the gas density.Left: The early slowing-down of the infall velocity at z ∼ 6 − 7 kpc coincides with a peak of the gas acceleration due to CRs with a particle momentum of ∼ 100 GeV c −1 .Right: The infalling gas has been converted to a coherent outflow up to z ≈ 22 kpc.The distribution of CRs extends to a similar height (z ∼ 25 kpc), which is a signature of advection dominated CR transport.This is also seen as a small enhancement in the density profile.The outflow front is driven by intermediate-to high-energy CRs, whereas the acceleration of the gas at lower heights is dominated by CRs with a particle momentum of p ∼ 10 GeV c −1 .
The time evolution is shown in Fig. 9 for the same four models as in Fig. 8.In the low-mass halo, CR transport is mainly dominated by advection for most of the time shown here.Only at the end, when no net outflows are driven, patches of outflowing gas are dominated by CR diffusion.The two intermediate-mass models indicate a transition from diffusion-dominated CR transport in the beginning to a more advective CR transport at later times.This corresponds to the initial infall, which is later converted to a coherent outflow.The outflows are strong and fast enough such that advective CR speeds are more important in comparison to CR diffusion.In contrast, the Milky Way analogue is dominated by diffusive speeds over the entire simulation time.
We further investigate the contribution of the individual CR particle energies to the dynamics.We do this by analysing vertical profiles in a cylinder of radius r = 15 kpc and a height of ±30 kpc that is centred on the galaxy.Figures 10  and 11 show the profiles for models "M 10 11 -spec-ΛA" and "M 10 12 -spec-ΛA", respectively.In each figure the left-hand panels show the profiles at t = 0.3 Gyr while the right-hand panels display the advancement of the outflow at t = 0.6 Gyr.The top two panels show profiles of the gas acceleration driven by CRs and the CR energy density, with the color of the lines indicating the CR momentum.The lower two panels show the vertical velocity and density, respectively.All panels show volume-weighted, azimuthally averaged quantities.
In the case of the 10 11 M⊙ halo (Fig. 10), the central region is dominated by low-energy CRs with p ≲ 30 GeV c −1 at both times.At the early stage, the dominant driver of the gaseous outflows are CRs with a momentum of p ∼ 100 GeV c −1 (via the energy-dependent CR acceleration |acr|).The peak acceleration coincides with both, an increase in velocity as well as a decrease in density.This indicates that this outflow has been launched from the centre of the galaxy.The acceleration profiles of the individual CR momenta are all spatially concentrated at around z = 6 kpc with a small spread of ∼ ±5 kpc.At later times, the inflow is reverted to an outflow, with the outflow front reaching a height of ∼ 22 kpc  10 but for the more massive galaxy (M 10 12 -spec-Λ A ) at t = 0.3 Gyr (left) and t = 0.6 Gyr (right).At both times we note a significantly larger spatial spread of the different CR energies in both the CR energy density as well as the CR-driven acceleration.This large spread is attributed to energy-dependent CR diffusion and highlights the diffusion dominated transport.The infall at t = 0.3 Gyr is significantly slowed down but not fully stopped by CRs with momenta of p ∼ 10 4 GeV c −1 for a height ≲ 10 kpc.The overall larger central densities result in more efficient CR cooling, in particular for low-energy CRs that cannot diffuse fast enough out of the dense gas.As a result, CRs with p ≲ 30 GeV c −1 neither dominate the energy nor the acceleration profile.at t = 0.6 Gyr.The CR energy density and the acceleration profiles have not advanced beyond the point of the actual outflow.This emphasizes the advection dominated transport of CRs, which implies that all CR momenta reach a similar distance from the centre.At first sight there is an inconsistency between the CR transport in model "M 10 11 -spec-ΛA" in Fig. 9, where we initially find predominantly CR diffusion, and Fig. 10, where advection clearly dominates for both times.This apparent discrepancy originates from the different volumes under consideration.The former investigates a larger cylindrical volume, which encompasses a significant fraction of the galactic disc.The latter analysis only focuses on the outflow cone with a strong direct CR driven impact.
The situation is different in the case of simulation "M 10 12 -spec-ΛA" (Fig. 11).At the early time (left-hand panels), the slowing down of the infall at z ∼ 10 kpc is mainly caused by CRs with high momenta (p ≳ 10 4 GeV c −1 ).We also note that there is a larger spread in the different CR momenta, corresponding to energy-dependent diffusion, which is typical for a diffusion dominated system.At later times, the infall is however not halted and the diffusive spread of CRs of different momenta is more pronounced.
By comparing the early times for both halo masses (Figs. 10 and 11) we note very similar maximum values of the magnitude of the outward pointing acceleration of CRs, |acr| ∼ 2 − 4 km s −1 Myr −2 .We therefore attribute the continued infall in the more massive galaxy mainly to the stronger gravitational attraction of the halo.In the shallower potential of the 10 11 M⊙ model, the same acceleration can revert the inflow to a net outflow.

TIME AVERAGED VERTICAL ACCELERATION
We further investigate the time averaged (0.5−2 Gyr) vertical accelerations for the individual components that are responsible for the dynamics.Again, we focus on the spectral models for each halo mass.We would like to stress that in order to identify the long-term evolution of the galaxy and a possible continuous wind, it is important to average the accelerations over a significant fraction of the simulation time.A single acceleration profile at one point in time could also indicate a short-term fountain flow effect rather than a long-lasting outflow.
We compute the outward pointing volume-weighted timeaveraged accelerations over 1.5 Gyr using 15 snapshots (Nsnap = 15).For the CR, and thermal acceleration we compute For the magnetic field, we need to take both magnetic pressure and tension into account, where the z-component [•]z is used for the plots.The volume average ⟨•⟩V is taken over segments of cylinder with height ∆z = 1 kpc and radius r, centred on the galaxy.First, we quantify the similarity of the CR accelerations.We show the time averaged outward pointing total CR acceleration in Fig. 12 for all halo masses for the spectral CR model.The halo with 10 12 M⊙ develops a five times stronger vertical CR acceleration compared to the 10 11 M⊙ halo, whereas the peak and average star formation rates are 13 and 9 times as high.The acceleration of the 3 × 10 11 M⊙ halo is twice as strong in comparison to the 10 11 M⊙ halo, which is similar to the scaling in star formation, where the peak and mean star formation rates are a factor of 2.8 higher in the larger halo.A much stronger difference is seen when comparing the 10 10 M⊙ halo with the 10 11 M⊙ counterpart.Here acr is lower by ∼ 30 per cent.The star formation rate is, however, more than an order of magnitude lower in the 10 10 M⊙ halo.
We further distinguish the time averaged accelerations between the different components in a cylinder with a radius of 3 kpc in Fig. 13.From top to bottom we increase the halo mass.In the left-hand panels we investigate the averaged accelerations from 0.5 − 1.25 Gyr while the right-hand counterpart shows the same quantity at a later time (1.25 − 2 Gyr).Shown are the total outward pointing acceleration (thermal, magnetic plus CR; dashed lines), the individual accelerations by CRs (solid lines), thermal pressure (dashed-dotted lines), and magnetic fields (dots).In addition, we plot the negative of the inward pointing gravitational attraction due to gas, stars and dark matter (dotted lines).In the case of the lowmass halo, the difference between the two averaging times does not significantly differ.Over the entire simulation time and at basically all heights the CR acceleration alone can compensate the gravitational attraction.The contribution of the other components (thermal and magnetic) is relatively small.The outflows driven from the 10 10 M⊙ halos are therefore basically purely CR driven, which is in line with the advection-dominated transport, see Section 7. The more massive models with 10 11 and 3 × 10 11 M⊙ show a different temporal evolution.Whereas in the beginning the CRs are the main driver of the outflow at all heights, at the later stage the net acceleration of the outflow is a combination of all energy components.The most massive model with 10 12 M⊙ clearly indicates the dominance of the gravitational attraction over the different outward pointing accelerations.The temporal differences are minor with only small spatial regions showing net outward pointing accelerations that exceed the gravitational attraction.The system is therefore dominated by locally and temporally driven fountain flows.The strong gravitational attraction furthermore forces the CR transport to be dominated by diffusion rather than advection.

CGM PROPERTIES
Finally, we investigate the CGM properties of our simulated galaxies with a focus on the thermal properties and the ratio of CR-to-thermal pressure.For this comparison, we also run galaxy models without CR feedback.Figures 14 and 15 illustrate the thermal properties for simulations with halo masses of 10 11 M⊙ and 10 12 M⊙, respectively, at t = 2 Gyr.From left to right we show models with different physics parameters, the panels from top to bottom depict edge-on views of the gas density, the temperature and the CR-to-thermal pressure ratio, Xcr, as well as a face-on view of Xcr.For the M = 10 11 M⊙ halo, the largest systematic difference is that the run without CRs shows a much smaller vertical extent of the disc in comparison to the runs including CRs.
The differences in the gas temperature map are even more apparent.Whereas the purely thermal simulation develops volume filling hot gas with T ∼ 10 6 K, the inclusion of CRs allows them to spread diffusively in the CGM, thereby providing pressure support and enabling the thermal gas to become thermally unstable and to cool down to temperatures of T ∼ 10 5 K.A noticeable exception are the outflow cones emerging from the centre of the galaxy.Here, the gas temperature exceeds 10 6 K for all CR models except the grey model with the highest diffusion coefficient (D = 3 × 10 29 cm 2 s −1 ).The difference in the CGM temperature between the grey and spectral models ("M 10 11 -grey4 × 10 28 -ΛA" vs. "M 10 11 -spec-ΛA") is marginal.Switching off Alfvén cooling increases the CR content in the outflow cone and reduces the gas temperature.In the grey model with the highest diffusion coefficient (D = 3 × 10 29 cm 2 s −1 ), the outflow cone is much less pronounced and the outflow region is colder.We also note morphological differences in the gas distribution.We attribute the much less pronounced outflow cone to the very fast diffusion, which in this case is larger than the energy weighted effective diffusion of the spectral model (cf.Fig. 7), and the resulting stronger CR-driven perturbations, which destroy the cone.
The differences between the CR models are even stronger in the ratio Xcr.Here, the CRs can quickly occupy the majority of the volume, which is dominated by CR pressure.In the grey CR model with D = 4 × 10 28 cm 2 s −1 , the outflow cone is marginally dominated by thermal pressure.Contrary, in the spectral model, the thermal pressure in the outflow cone is an order of magnitude larger than the CR pressure.This difference is a consequence of the energy dependent diffusion.
Effectively, more CR energy can escape via fast diffusion, which lowers the CR content in the outflow while the ther- mal structure of the outflow cone is comparable.We note that the advection dominated CR transport discussed in Section 7 focuses on early times and evolves to a balanced transport situation after t ≳ 1 Gyr (Fig 9).When Alfvén cooling is switched off, the less efficient cooling of CRs results in equal pressures of CRs and thermal gas in the outflow cone.
For the Milky Way-mass galaxy model (Fig. 15), the situation is different because no significant net outflows are launched.The purely thermal simulation forms a very thin disc with no perturbations driven into the halo.This implies a hot CGM with temperatures exceeding 10 6 K that has been heated by the outwards propagating accretion shock.Including CRs allows for a colder environment at the disc-halo interface.This is due to CR-driven fountain flows from the galactic disc, which entrain the initially hot CGM and mix it with the colder ISM, thereby enabling thermal instability.As a result, the CRs are providing the dominant pressure support at the halo-disc interface.In this case, the CR transport speed is of importance.The grey run with a low diffusivity ("M 10 12 -grey4 × 10 28 -ΛA") drives weaker fountain flows that reach a lower altitude above the disc in comparison to the grey simulation with larger diffusion coefficients (M 10 12 -grey10 29 -ΛA and M 10 12 -grey3 × 10 29 -ΛA).The fast diffusion of the high-energy CRs in the spectral model (M 10 12 -spec-ΛA) is comparable to the grey model with D = 10 29 cm 2 s −1 , which is in agreement with the energy  (right), measured in a cylinder of radius 25 kpc and height ±20 kpc that is centred on the galaxy.In all cases the purely thermal runs have the highest temperature in the inspected volume.For the two lowest mass halos, the temperatures in the CR runs do not differ significantly.For the two highest mass halos, the grey runs with large diffusion coefficient and the spectral model show noticeably lower temperatures.For Xcr there is a general trend for all halo masses: Xcr increases for the grey models with increasing diffusivity.The values for the spectral models depend on their effective diffusion coefficient.
weighted effective diffusion coefficient, see Fig. 7.The strong impact of faster diffusion on the thermal properties of the (inner) CGM is also seen in the face-on view of Xcr (bottom panels).Whereas the grey models with D = 4 × 10 28 and 10 29 cm 2 s −1 show a sharp edge at r ∼ 40 kpc beyond which no significant amount of CRs can be transported, the high-D grey and the spectral run fill a larger volume with CR pressure dominated gas (Xcr > 10).
We show the time averaged values for the temperature and Xcr over 1.5 Gyr (t = 1 − 2.5 Gyr) in Fig. 16.Both measures are volume weighted in a cylinder of radius r = 25 kpc and height z = ±20 kpc around the centre of the galaxy.We define ⟨Xcr⟩V = ⟨Pcr/P th ⟩V as the volume-weighted CRto-thermal pressure ratio.As illustrated in Figs. 14 and 15, the simulations without CRs form a hotter CGM.The gas temperature increases with halo mass, which is a direct consequence of a virialisation process of the CGM driven by the expanding accretion shock.In all but the Milky Way models the difference between the non-CR and CR runs is more than a factor of two in the average temperature.In case of the two low-mass halos, the differences between the different CR physics does not alter the temperatures significantly.This is in line with the temporally strongly advective motions and the resulting outflows, where the injected CRs and the cooler gas are equally well transported throughout the inspected volume.For the two more massive halos CR diffusion plays a more important role, which is directly reflected in ⟨T ⟩.The grey models with lower diffusivity have a weaker impact on the halo gas and lead to a higher gas temperature in comparison to the grey models with a larger diffusion coefficient and the spectral model.For ⟨Xcr⟩ we first note that in all halos, the CR pressure dominates over the ther-mal pressure (in the considered volume) with a smallest value of ⟨Xcr⟩ ∼ 5. We observe an increasing value of ⟨Xcr⟩ with increasing diffusivity for all halo masses, which is exactly reversed in comparison to the temperature, just as we expect if the halos (outside the fast outflow) are in approximate hydrostatic equilibrium.The lowest value of ⟨Xcr⟩ is measured for models with D = 10 28 cm 2 s −1 .The corresponding runs without Alfvén cooling leave more CR energy in the system, which increases ⟨Xcr⟩.

DISCUSSION
We find that between a halo mass of 3×10 11 and 10 12 M⊙ the gravitational attraction towards the centre of mass becomes too strong for the stellar feedback to drive outflows in our setting.This result is consistent with previous studies.Booth et al. (2013) model stellar feedback including CRs in both a dwarf as well as a Milky Way galaxy and find that only the former one is able to drive outflows with mass loading factors exceeding unity.A similar limiting mass for driving sustained outflows over longer times was found in CR models by Jacob et al. (2018) using an overall very similar setup like our grey models, except that they use a Hernquist (Hernquist 1990) instead of an NFW profile (Navarro et al. 1997) and isotropic CR diffusion.However, the details in the star formation rate and the mass loading factor differ.We find mass loading factors that are approximately η ∼ 0.5 − 2 for the 10 10 M⊙ halo and η ∼ 0.5 for the 10 10 M⊙ model, whereas Jacob et al. (2018) record η ∼ 20 − 40 and η ∼ 2 − 3 for the same halo masses mass.We note that the models that have investigated the halo mass dependence only use a sim-plified model for the stellar feedback.Small scale simulations of the multiphase ISM show that the details of early stellar feedback like radiation and stellar winds in combination with SN feedback and CRs can perceptibly change the effective outflows significantly (e.g.Gatto et al. 2017;Rathjen et al. 2021;Simpson et al. 2022;Kim et al. 2022).This is particularly important in low-mass galaxies, where stellar feedback has a strong effect.
Besides the ability to drive net outflows, CRs change the thermal properties in the ISM and the outflow.Girichidis et al. (2018) find smoother and colder gas in the presence of CRs compared to non-CR models using stratified boxes.This is consistent with full galactic models (Dashyan & Dubois 2020;Nuñez-Castiñeyra et al. 2022).Apart from the thermal structure, CRs leave an imprint on the angular momentum of the outflowing gas.Peschken et al. (2021Peschken et al. ( , 2022) ) find that there are two different outflows phases; the warm outflows with high angular momentum and the hot outflows with low angular momentum.
We chose the simplified advection-diffusion approximation for the CR transport, which does not explicitly model the interaction of CRs with Alfvén modes due to the streaming instability (Kulsrud & Pearce 1969;Wentzel 1974;Skilling 1971Skilling , 1975;;Shalaby et al. 2021).The amplification of Alfvén waves requires energy, which is provided by the CRs and causes them to lose energy.Since diffusion is energy conserving, this loss of CR energy is not included in the pure CR diffusion approach.The effect of CR losses during CR streaming can to first order be emulated in a modified diffusion approach, (Wiener et al. 2013a;Buck et al. 2020), which is also included in this study.This means that our effective transport speeds implicitly depend on the effective Alfvén wave cooling since it alters the CR energy gradients.Sharma et al. (2010) included CR streaming in a one-moment approach (see also Wiener et al. 2017;Holguin et al. 2019), in which the streaming speed needs to be regularised.However, this introduces numerical parameter that is connected to a critical CR gradient length and causes an unphysically fast transport of CRs if they injected into an already existing distribution of background CRs (Thomas & Pfrommer 2019).
A more fundamental approach uses two moments to compute the CR transport and the resulting energy losses.Jiang & Oh (2018); Chan et al. (2019); Armillotta et al. (2021Armillotta et al. ( , 2022) ) evolve the CR energy density and the momentum density with a diffusion coefficient that assumes a steady-state wave energy.More accurately, Thomas & Pfrommer (2019, 2022) solve for the forward and backward propagating resonant Alfvén wave energies, which allows to compute the CR diffusion coefficient more self-consistently.In these twomoment approaches the effective CR transport speeds directly depend on the efficiency of wave damping and the coupling of the CRs to the ISM.While this approach is more self-consistent, it requires an estimate of the degree of ionization of the background plasma since the efficiencies of ionneutral and non-linear Landau damping of Alfvén waves are strong functions of the ionization degree.The resulting effective transport speed therefore sensitively depends on an accurate modeling of the different ISM phases.We note that the numerical setup of an isolated cooling halo is ideal for systematically investigating the internal processes in the forming galaxy.However, none of the more complex environmental conditions expected in cosmological simulations are accounted for, such as the long-term evolution including accretion of gas onto the halo, development of turbulent motions and the accompanied amplification of magnetic fields in the halo (e.g.Pakmor et al. 2020).This emphasizes the importance of applying the spectral CR model to cosmological (zoom) setups.

CONCLUSIONS
We have performed simulations of isolated rotating cooling halos of different masses from dwarfs to Milky Way analogues with four masses (M = 10 10 , 10 11 , 3 × 10 11 , and 10 12 M⊙).The forming galaxies use an effective model for the ISM including a stochastic star formation and SN feedback prescription, in which we only inject the CR energy for each SN directly.We inject 10 per cent of the SN energy in CRs, which we model as a dynamically coupled fluid.We distinguish between a grey CR approximation, in which only the total CR energy is evolved and a spectrally resolved counterpart, where we resolve the CR momenta ranging from the subrelativistic (pmin = 0.1 GeV c −1 ) up to the ultra-relativistic regime (pmax = 100 TeV c −1 ).We use the advection-diffusion approximation for the CR transport with anisotropic energy dependent diffusion along the magnetic field lines.In order to account for CR losses during the transport we also apply a simplified model for CR Alfvén wave cooling.We perform 32 simulations in total, which allow us to draw the following conclusions.
• We find CR driven outflows in all but the most massive model, i.e., for all halos below the Milky Way analogue.The overall mass loading factors reduce with increasing halo mass due to the stronger gravitational attraction.In the Milky Way model (M = 10 12 M⊙) the strong initial infall of halo gas slows down over time but cannot be reverted to a net outflow.The time averaged acceleration profiles reveal that in the lowmass model CRs alone are able to drive the outflows.In the intermediate mass halos the outflows are dominantly but not exclusively driven by CRs.
• We find a clear distinction between advective and diffusive CR transport.In the galaxy with the lowest halo mass, the low gravitational impact allows the CRs to have a strong dynamical effect.Strong outflows with large mass loading factors are driven shortly after the onset of star formation.As a result, the transport of CRs is mainly advective in this model.The three more massive halos experience a stronger initial infall of gas before feedback can slow down the infall.During this initial phase the three models show a more diffusive CR transport.For the two intermediate mass halos (M = 10 11 M⊙ and M = 3 × 10 11 M⊙) the onset of a delayed outflow can transform the overall diffusive to a globally advective CR transport.The outflow cone is dominated by advective transport from the very beginning.In the Milky Way model, the transport is dominated by diffusion over the entire simulated time.We clearly identify typical diffusive features where higher energy CRs reach larger distances from the injection sites despite a net infall of gas.
• The effective diffusion coefficients in the spectral models vary by two orders of magnitude in space and over time.In the injection regions close to the centres of the galaxies, where most stars form, the spectra are close to the injection spectra.Since we inject CRs with a power-law with f (p) ∝ p −4.2 the spectra are dominated by low energy CRs.The corresponding energy weighted diffusion coefficients are therefore smaller close to the injection sites (D(p) ∝ p 0.3 ) with values of ⟨D⟩e cr ∼ 10 28 cm 2 s −1 .At larger distances from the star forming centres, the effective diffusion coefficients depend on the dominant CR transport mode.For the models with M = 10 10 M⊙ and M = 10 11 M⊙ strong outflows lead to a predominantly advective CR transport, which keeps the shape of the spectrum and results in similarly low values of ⟨D⟩e cr ∼ 10 28 cm 2 s −1 in the outflow region.For the Milky Way model, the predominantly diffusive transport causes the spectra in outer the regions to be dominated by CRs with higher momenta.The resulting values for the effective diffusivity are then ⟨D⟩e cr ∼ 2 − 3 × 10 29 cm 2 s −1 .
• The scaling of the CR energy density with the gas density closely follows an adiabatic scaling of ecr ∝ ρ 4/3 for densities ρ ≲ 10 −26 g cm −3 .Above this density, non-adiabatic effects such as CR cooling and CR diffusion out of the dense gas flatten the scaling.Hadronic cooling scales with Λ hadr ∝ ρ 7/3 for the density regime of adiabatic CRs.Alfvén cooling scales with the density as predicted by simple scaling estimates with ΛA ∝ ρ 11/6 .Except for the lowest densities (ρ ≲ 10 −28 g cm −3 ) the hadronic cooling dominates over the Alfvénic counterpart, which explains the small differences in star formation and outflow rate between simulations with and without Alfvén cooling.
• Both star formation and outflow rates do not vary perceptibly between the different CR models.In particular, spectrally resolved CRs do not change the peak of the star formation rate or the overall temporal evolution 6 .Minor differences are present for the lowest mass halo (M = 10 10 M⊙) due to the overall large dynamical impact of the CRs.This demonstrates the need to account for spectral CR transport in order to provide a reliable description of the CGM as well as CR-induced non-thermal radio and gammaray emission (Werhahn et al. 2023). .Effective adiabatic index γcr (left), the total CR pressure Pcr (middle) and the pressure ratio Pcr/P cr,simple (right), where P cr,simple is computed using γ = 4/3.Here we show simulation "M 10 12 -spec-Λ A " at t = 2 Gyr.The other simulations show comparable or even smaller differences between the spectrally computed pressure and the simplified counterpart.

APPENDIX B: EMULATED CR STREAMING LOSSES
Modelling the losses due to the streaming instability should depend on the individual CR energies.In particular, the high energy CRs with momenta p ≳ 1 TeV c −1 are likely not affected by the streaming losses.Our emulated streaming losses reduce the net CR energy and we simply rescale the spectrum to account for the total desired loss.This means that the effective losses in each spectral bin depend on the spectral shape, because the peaks in the spectra suffer a larger net loss compared to bins which only carry a small amount of CR energy.For most cases, the spectra peak in the range between 1 − 10 GeV c −1 .The critical spectral bins (bin numbers 8-12, p > 1 TeV c −1 ), only contain a small amount of total CR energy in each cell, so the error in applying a simplified cooling are expected to be small.Figure B1 shows the fraction of the energy in each spectral bin for simulation "M 10 11 -spec-ΛA" at t = 2 Gyr.Shown are the radially averaged spectra (lefthand panel) as well as the vertically averaged counterparts (right-hand panel) for different concentric rings and heights, respectively.

APPENDIX C: MAGNETIC FIELD STRUCTURE
We illustrate the magnetic morphology via magnetic field lines for all spectral models in Fig. C1.Depicted is the gas density color coded with overplotted streamlines.While being predominantly toroidal, the fields in the disc show more small scale features in comparison to the CGM and the outflow region, except for the smallest halo mass.The variations between different CR models for each halo mass are subdominant.

−
Figure1.Distribution of the CR energy density, the hadronic and the Alfvén cooling rate as a function of density in the simulation "M 10 11 -spec-Λ A " averaged over t = 1 − 2 Gyr.Colour coded is the distribution in these variables, weighted with the fraction of CR energy density.Overlaid are the median (dashed lines) as well as the analytical scaling that assumes adiabatic CRs only (solid lines), which is only valid at low densities.Hadronic cooling dominates over Alfvén cooling above a density of ρ ∼ 10 −28 g cm −3 .This explains the small differences in the galaxy properties between the simulations with and without Alfvén cooling.

Figure 3 .
Figure3.Kennicutt-Schmidt relation for all simulations, colour coded by their halo masses (see legend).The differences between the models at the identical halo mass is negligible.The different halo masses mainly occupy different ranges in the gas column density.Below a surface density of Σgas ≈ 10 M ⊙ pc −2 the star formation rate reduces compared to the classical relation, which is a well-known observed feature (e.g.,Bigiel et al. 2008;Leroy et al. 2008) . The mass flux is measured at a height of |z| = 10 kpc through a circle with radius r = 15 kpc.The initial collapse of the gas is imprinted in the early

Figure 5 .
Figure5.Time-averaged (0.5 − 2.5 Gyr) mass loading factor as a function of total mass.The outflow is measured as in Fig.4and the mass loading factor is computed for outflow and star formation rates at the same time.All but the highest mass galaxy develop effective outflows with positive mass loading factors.The differences between the CR models is stronger for lower masses.The largest mass loading factors are reached by the grey simulations with D = 4 × 10 28 and 10 29 cm 2 s −1 .The temporal scatter is indicated by the errorbars marking the 25 and 75 percentile.

MFigure 6 .
Figure6.Onset of the outflow.The four left-hand columns depict the grey models with increasing diffusion coefficient "M 10 12 -grey10 28 -Λ A ", "M 10 12 -grey4 × 10 28 -Λ A ", "M 10 12 -grey10 29 -Λ A ", "M 10 12 -grey3 × 10 29 -Λ A " while the right-hand column shows the spectral model "M 10 12 -spec-Λ A ".The top two rows show the edge-on velocity and density structure while the two bottom rows display the face-on density and total CR energy density, respectively.With increasing diffusion coefficient, the outflow front expands faster.At this early stage the spectral model behaves similarly as the grey simulation with D = 10 29 cm 2 s −1 .The main structure of the disc is very similar for most cases; only the grey model with D = 3 × 10 29 cm 2 s −1 forms a noticeably smaller disk.

MM
Figure7.Effective CR energy weighted diffusion coefficient at t = 0.2 Gyr (top) and t = 1 Gyr (bottom) for models with D 0 = 10 28 cm 2 s −1 .From left to right we show cuts through the middle of the box for halo masses of 10 10 , 10 11 , 3 × 10 11 , and 10 12 M ⊙ .We exclude regions that are not yet populated by CRs (white areas).The circles indicate the radius of 75 per cent enclosed stellar mass.At early times we find small effective diffusivities (⟨D⟩e cr ∼ 10 28 cm 2 s −1 ) close to the centre, where new CRs are injected.At larger distances from the centre the spectra are more dominated by high energy CRs, which increases ⟨D⟩e cr to values of ∼ 3 × 10 29 cm 2 s −1 .At late times we note a strong difference between the regions that are dominated by advective transport with low-energy CRs dominating the spectrum and regions of diffusive transport with higher diffusivities.

MFigure 8 .
Figure 8.Comparison of the advective and diffusive CR speeds in the spectral CR models at t = 0.5 Gyr.From left to right we increase the halo mass.From top to bottom we show edge-on maps of the gas density, the radially outward pointing advection velocity (vr), the diffusion speed (v diff ) as well as their ratio vr/v diff .The low mass halo develops a smooth coherent advective outflow.The two medium mass systems show patches of strong advective motions.These features are very weak in the infall dominated high-mass halo.The ratio of the two speeds indicates that advective velocities are more dominant in low-mass systems, whereas high-mass systems are dominated by CR diffusion, see also Fig. 9.

MFigure 9 .
Figure9.Volume averaged ratio ⟨vr/v diff ⟩ V within a cylinder of radius r = 25 kpc and height ±12.5 kpc centred on the galactic centre.From low to high masses there is a transition from advectively to diffusively dominated CR transport.The two intermediate mass halos start with diffusive CR transport before strong continuous CR-driven outflows develop so that CRs are primarily transported via advection.
Figure10.Vertical profiles averaged over a cylinder with radius r = 15 kpc and height ±30 kpc positioned on the centre of the galaxy at t = 0.3 Gyr (left) and t = 0.6 Gyr (right) for simulation "M 10 11 -spec-Λ A ". From top to bottom we show the acceleration of the gas owing to CRs, their energy density, the gas velocity as well as the gas density.Left: The early slowing-down of the infall velocity at z ∼ 6 − 7 kpc coincides with a peak of the gas acceleration due to CRs with a particle momentum of ∼ 100 GeV c −1 .Right: The infalling gas has been converted to a coherent outflow up to z ≈ 22 kpc.The distribution of CRs extends to a similar height (z ∼ 25 kpc), which is a signature of advection dominated CR transport.This is also seen as a small enhancement in the density profile.The outflow front is driven by intermediate-to high-energy CRs, whereas the acceleration of the gas at lower heights is dominated by CRs with a particle momentum of p ∼ 10 GeV c −1 .

Figure 11 .
Figure11.Same as Fig.10but for the more massive galaxy (M 10 12 -spec-Λ A ) at t = 0.3 Gyr (left) and t = 0.6 Gyr (right).At both times we note a significantly larger spatial spread of the different CR energies in both the CR energy density as well as the CR-driven acceleration.This large spread is attributed to energy-dependent CR diffusion and highlights the diffusion dominated transport.The infall at t = 0.3 Gyr is significantly slowed down but not fully stopped by CRs with momenta of p ∼ 10 4 GeV c −1 for a height ≲ 10 kpc.The overall larger central densities result in more efficient CR cooling, in particular for low-energy CRs that cannot diffuse fast enough out of the dense gas.As a result, CRs with p ≲ 30 GeV c −1 neither dominate the energy nor the acceleration profile.

Figure 12 .
Figure12.Comparison of time averaged (0.5 − 2 Gyr) vertical profiles of the CR-driven acceleration for all halo masses.We use the spectral CR model for each halo and the profiles are averaged within a cylinder of radius r = 15 kpc centred on the galaxy.The differences between the halo masses are smaller than the differences in the star formation and CR injection rate, see text.

Figure 13 .MFigure 14 .
Figure13.Time averaged vertical acceleration profiles in a cylinder with radius r = 3 kpc for the spectral CR model and all halo masses from top to bottom.The left-hand panels show the average from 0.5 − 1.25 Gyr, the right-hand ones for 1.25 − 2 Gyr.The low mass halo with 10 10 M ⊙ is dominated by CR acceleration, which explains the strong outflows with high mass loading (see Fig.5).The two intermediate mass systems show similarities in their profiles with CRs being dominant in the beginning and supportive but not exclusively responsible for the outflow at later stages.In the most massive halo the outward pointing forces cannot compensate the gravitational attraction, i.e., no wind is launched from this galaxy.

MFigure 15 .
Figure15.Same as Fig.14for models with a halo mass of 10 12 M ⊙ at t = 2.0 Gyr.There are no significant net outflows launched for this halo mass.Nonetheless small CR-driven fountain flows cause the CGM to cool to lower temperatures in models including CRs.As a result, the cold regions are characterised by a dominating CR pressure rather than by thermal pressure.

Figure 16 .
Figure16.Time averaged (t = 1 − Gyr) volume weighted values of gas temperature, T (left), and CR-to-thermal pressure ratio, Xcr (right), measured in a cylinder of radius 25 kpc and height ±20 kpc that is centred on the galaxy.In all cases the purely thermal runs have the highest temperature in the inspected volume.For the two lowest mass halos, the temperatures in the CR runs do not differ significantly.For the two highest mass halos, the grey runs with large diffusion coefficient and the spectral model show noticeably lower temperatures.For Xcr there is a general trend for all halo masses: Xcr increases for the grey models with increasing diffusivity.The values for the spectral models depend on their effective diffusion coefficient.
Figure A1.Effective adiabatic index γcr (left), the total CR pressure Pcr (middle) and the pressure ratio Pcr/P cr,simple (right), where P cr,simple is computed using γ = 4/3.Here we show simulation "M 10 12 -spec-Λ A " at t = 2 Gyr.The other simulations show comparable or even smaller differences between the spectrally computed pressure and the simplified counterpart.

Table 2 .
Spectral bins and the corresponding diffusion coefficients.Listed are the bin number, the left-hand momentum boundary, the central momentum as well as the corresponding diffusion coefficient for three different scalings of D with the momentum.